Re:

[Konstantinos Poulios]

Dear Anne-Cecile

I had responded to your original mail (25 oct.) with the following answer, but it probably never arrived to you due to email filter not allowing python files. Trying again.

A.
Yes the error message makes sense, your mfp mesh_fem is reduced because of the line

mfp = gf.MeshFem("partial", mfp_, keptdofs)

You should instead use:

md.set_variable("p", mfp.reduce_vector(md.interpolation("p0*min((X(2)-0)/10,1)", mfp_)))


B.
I am not sure what you mean with implementing Eq. (12)

This equation looks like part of the derivation of the analytical solution that you do not need for your numerical solution. You actually need to compare your numerical solution to this result, if that's what you mean, you can compute this quantity in the whole domain and export it with something like
mfout.export_to_vtu( ......., mfout, md.interpolation("1/a*Grad_u(2,2)-alpha*p",mfout), "Sigma_yy_from_eq12")

In general your benchmark seems to be 1D, and you are solving it numerically in 2D. So you just have to verify that the computed fields are constant along the x-direction.

What you actually miss in your script is:
#Equation 3 sigma_y =-p0 loading #missing
md.add_linear_term(mim4, 'p0*u(2)', T_RG )

which will impose the sigma_y=-p0 condition at the top (this is quite standard in continuum mechanics, mathematically proven by use of the divergence theorem)

C.
dp/dn=0 should be true except for the top boundary where you impose p=0


D.
I have also found some mistakes in your script:
- For the definition of the boundary regions, the direction vectors should point outwards.
- You defined the terms "multL*u(1)", "multR*u(1)" and "multB*u(2)" in the whole domain instead of the boundaries. Not sure why GetFEM did not give you an error there.

After these fixes, the attached version og the script seems to run correctly. I have changed some input parameters compared to your version and added some more output. Use meld (https://meldmerge.org/) to compare with your version.

Best regards
Konstantinos

On Mon, Nov 8, 2021 at 7:52 PM Lesage,Anne Cecile J <AJLesage@mdanderson.org> wrote:

Dear getfem users

 

I am trying to implement the Bio equation for brain shift modeling during surgery

I am searching to solve a benchmark in 2d and 3d with a 1d analytical solution see attached file

I have a message error in python when I try to interpolate my initial pressure field see attached picture

The error is in the line

md.set_variable("p",md.interpolation("p0*min((X(2)-1)/10,1)",mfp)) 

In this first example my mesh is triangular 2d

 

Thank you

Anne-Cecile

 

#!/usr/bin/env python3

# -*- coding: UTF8 -*-

#

 

import sys

sys.path.append('/usr/local/lib/python3.8/site-packages/getfem')

 

import getfem as gf

 

import numpy as np

gf.util_trace_level(1)

 

# Input data

G = 1e7         # [N/m^2]

nu = 0.3

k = 1e-13       # [m^3 s/kg]

dt = 1e3        # [s]

steps = 1000

alpha = 1.0        # ratio of fluid volume extracted to volume change of the tissue under compression

invs = 0

 

p0 = 1e3        # [N/m^2] normal traction at the top

 

press_fem_order = 1      # pressure finite element order

disp_fem_order = 2       # displacements finite element order

mult_fem_order = 2       # displacement multipliers finite element order

 

 

 

print('Read Mesh GiD')

 

#gf.util('trace level', 2)   # No trace for mesh generation

#mesh = gf.Mesh('generate', mo, h, 2)

 

mesh=gf.Mesh('import','gid','mesh2d_h0_2.msh')

 

 

   

# mesh generation

#NX = 2*((NX+1)//2)

#NY = 2*((NY+1)//2)

#XX = np.linspace(-LX/2, LX/2, NX+1)

#YY = np.linspace(0, LY, NY+1)

#mesh = gf.Mesh("cartesian", XX, YY)   

 

export_mesh = True        # Draw the mesh after mesh generation or not

 

# print mesh vtk format#

if (export_mesh):

    mesh.export_to_vtk('mesh.vtk');

    print('\nYou can view the mesh for instance with');

    print('mayavi2 -d mesh.vtk -f ExtractEdges -m Surface \n');

 

# mesh regions

#L_RG = 6

#B_RG = 7

#R_RG = 8

#T_RG = 9

#m.set_region(L_RG, m.outer_faces_with_direction([-1,0],0.01)

#m.set_region(B_RG, m.outer_faces_with_direction([0,-1],0.01)

#m.set_region(R_RG, m.outer_faces_with_direction([1,0],0.01)

#m.set_region(T_RG, m.outer_faces_with_direction([0,1],0.01)

 

#

# Boundary selection

#

# mesh regions

L_RG = 6

B_RG = 7

R_RG = 8

T_RG = 9

fb1 = mesh.outer_faces_with_direction([ 1., 0.], 0.01) # Right boundary

fb2 = mesh.outer_faces_with_direction([-1., 0.], 0.01) # Left boundary

fb3 = mesh.outer_faces_with_direction([0., -1.], 0.01) # Top boundary

fb4 = mesh.outer_faces_with_direction([0., 1.], 0.01) # Bottom boundary

 

RIGHT_BOUND=8; LEFT_BOUND=6; TOP_BOUND=9; BOTTOM_BOUND=7;

 

# Define mesh region

mesh.set_region( R_RG, fb1)

mesh.set_region( L_RG, fb2)

mesh.set_region( T_RG, fb3)

mesh.set_region( B_RG, fb4)

 

#mesh.region_merge(LATERAL_BOUND, RIGHT_BOUND)

 

 

# FEM

mfp_ = gf.MeshFem(mesh, 1)

mfp_.set_classical_fem(press_fem_order)

keptdofs = np.arange(mfp_.nbdof())

keptdofs = np.setdiff1d(keptdofs, mfp_.basic_dof_on_region(T_RG))

mfp = gf.MeshFem("partial", mfp_, keptdofs)

 

mfu = gf.MeshFem(mesh, 2)

mfu.set_classical_fem(disp_fem_order)

 

mfmult = gf.MeshFem(mesh, 1)

mfmult.set_classical_fem(mult_fem_order)

 

mfout = gf.MeshFem(mesh)

mfout.set_classical_discontinuous_fem(2)

 

# Integration methods

mim4 = gf.MeshIm(mesh, 3)

mim9 = gf.MeshIm(mesh, 5)

#mimd4 = gf.MeshImData(mim4, -1) # use these for saving data on integration points

#mimd9 = gf.MeshImData(mim9, -1) #

 

# Model

md = gf.Model("real")

md.add_fem_variable("u", mfu)      # displacements field

md.add_fem_variable("p", mfp)      # hydrostatic pressure field

 

md.add_filtered_fem_variable("multL", mfmult, L_RG)

md.add_filtered_fem_variable("multR", mfmult, R_RG)

md.add_filtered_fem_variable("multB", mfmult, B_RG)

 

md.add_fem_data("u_prev", mfu)

md.add_fem_data("p_prev", mfp)

 

md.add_initialized_data("G", G)

md.add_initialized_data("nu", nu)

md.add_initialized_data("k", k)

md.add_initialized_data("alpha", alpha)

md.add_initialized_data("invs", invs)

md.add_initialized_data("dt", dt)

md.add_initialized_data("p0", p0)

 

#Equation 1

md.add_linear_term(mim9, 'G*Grad(u):Grad(Test_u)+G/(1-2*nu)*Div(u)*Div(Test_u)+alpha*Grad(p).Test_u')

 

#Equation 2

md.add_linear_term(mim4, 'alpha/dt*(Div(u)-Div(u_prev))*Test_p+k*Grad(p).Grad(Test_p)')

 

#Equation 3 sigma_y =-p0 loading #missing

 

###############################

# initial and boundary conditions

#######################################

# md.add_linear_term(mim9, "p0*Test_u(2)", T_RG) ## sigma_n = pO???

md.add_linear_term(mim9, "multL*u(1)")  # u_n =0

md.add_linear_term(mim9, "multR*u(1)")  # u_n =0

md.add_linear_term(mim9, "multB*u(2)")  # u_n =0

 

# Neumann homogeneous on pressure, dp/dn=0 on L,R and B

 

# p=0 at T_RG is already enforced in the way mfp is defined from mfp_

md.set_variable("p",md.interpolation("p0*min((X(2)-1)/10,1)",mfp))  #  Initial condition p=p0

 

print('elements number=%d, points number=%d' % \

      (mesh.nbcvs(),mesh.nbpts()))

 

for step in range(steps):

   nit, conv =  md.solve("noisy", "lsolver", "mumps", "max_iter", 20, "max_res", 1e-9,

                         "lsearch", "simplest", "alpha max ratio", 1e9, "alpha min", 1.,

                         "alpha mult", 0.1, "alpha threshold res", 1e9)

   if(step%10==0):     

     print('\n print results every 10 time step');               

     mfout.export_to_vtu("Biot_benchmark_%i.vtu" % step,

                       mfu, md.variable("u"), "Displacements",

                       mfp, md.variable("p"), "Pressure")

     md.set_variable("u_prev", md.variable("u"))

     md.set_variable("p_prev", md.variable("p"))

The information contained in this e-mail message may be privileged, confidential, and/or protected from disclosure. This e-mail message may contain protected health information (PHI); dissemination of PHI should comply with applicable federal and state laws. If you are not the intended recipient, or an authorized representative of the intended recipient, any further review, disclosure, use, dissemination, distribution, or copying of this message or any attachment (or the information contained therein) is strictly prohibited. If you think that you have received this e-mail message in error, please notify the sender by return e-mail and delete all references to it and its contents from your systems.

attachments.zip
Description: Zip archive