[no subject]

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"))

From:	Lesage,Anne Cecile J
Subject:
Date:	Mon, 8 Nov 2021 18:52:14 +0000