Dear all,
I am using getfem with python interface (svn revision 4840). In my program the
level set changes at each iteration. Aftrewards i have to update the
"mesh_im_level_set " for different level set, but unfortunately the program
crashs and gives error like "Segmentation fault" (see below). For more details,
I have a sample of the program enclosed. If someone can help i will be grateful.
Thank you in advance.
Yao Koutsawa---------------------------------------->Adapting the mesh
Integration methods adapted
Trace 2 in getfem_models.cc, line 6448: Stiffness matrix assembly for isotropic
linearized elasticity
Trace 2 in getfem_models.cc, line 6448: Stiffness matrix assembly for isotropic
linearized elasticity
Trace 2 in getfem_models.cc, line 3938: Mass term assembly for Dirichlet
condition
Trace 2 in getfem_models.cc, line 3626: Source term assembly
iter 1
Integration methods adapted
Trace 2 in getfem_models.cc, line 6448: Stiffness matrix assembly for isotropic
linearized elasticity
Segmentation fault<----------------------------------------
import getfem as gf
import numpy as np
from numpy import sin, cos, pi
NY = 50
ls_deg = 1
deg = 1
# Mesh definition, simplices or square or ....
mesh = gf.Mesh("regular simplices", np.linspace(-1.0,1.0, 2*NY+1),
np.linspace(-0.5,0.5, NY+1))
# boundary selection
flst = mesh.outer_faces()
fnor = mesh.normal_of_faces(flst)
tleft = abs(fnor[1,:]+1) < 1e-14
ttop = abs(fnor[0,:]-1) < 1e-14
fleft = np.compress(tleft, flst, axis=1)
ftop = np.compress(ttop, flst, axis=1)
# mark it as boundary
DIRICHLET_BOUNDARY_NUM1 = 1
DIRICHLET_BOUNDARY_NUM2 = 2
NEUMANN_BOUNDARY_NUM = 3
mesh.set_region(DIRICHLET_BOUNDARY_NUM1, fleft)
mesh.set_region(DIRICHLET_BOUNDARY_NUM2, ftop)
# define level set on the defined mesh
ls = gf.LevelSet(mesh, ls_deg)
mf_ls = ls.mf()
ULS = mf_ls.eval('(-0.6+sin(pi*4*x)*cos(pi*4*y))', globals(),locals())
ls.set_values(ULS)
# the mesh with LS to reprensent the cut
mesh_ls = gf.MeshLevelSet(mesh)
mesh_ls.add(ls)
print('Adapting the mesh')
mesh_ls.adapt()
mim_in = gf.MeshIm('levelset', mesh_ls,'inside', gf.Integ('IM_STRUCTURED_COMPOSITE(IM_TRIANGLE(6),3)'))
mim_in.set_integ(4)
mim_out = gf.MeshIm('levelset', mesh_ls,'outside', gf.Integ('IM_STRUCTURED_COMPOSITE(IM_TRIANGLE(6),3)'))
mim_out.set_integ(4)
mf_basic = gf.MeshFem(mesh)
mf_basic.set_fem(gf.Fem("FEM_PK(%d,%d)"%(2,deg)))
mf = gf.MeshFem('levelset', mesh_ls, mf_basic)
mf.set_qdim(2)
print('Integration methods adapted')
# Model
md = gf.Model("real")
md.add_fem_variable("u", mf)
md.add_initialized_data('Force', [0.0, 1.0])
md.add_initialized_data('lambda1', 1)
md.add_initialized_data('mu1', 1)
md.add_initialized_data('lambda2', 2)
md.add_initialized_data('mu2', 2)
md.add_isotropic_linearized_elasticity_brick(mim_in, "u", "lambda1", "mu1")
md.add_isotropic_linearized_elasticity_brick(mim_out, "u", "lambda2", "mu2")
md.add_Dirichlet_condition_with_penalization(mim_in, 'u', 1.E+9, DIRICHLET_BOUNDARY_NUM1)
md.add_source_term_brick(mim_in, 'u', 'Force', DIRICHLET_BOUNDARY_NUM2)
for i in xrange(10):
if(i > 0):
ULS = ULS + 0.0025*i
ls.set_values(ULS)
mesh_ls.adapt()
mim_in.adapt()
mim_out.adapt()
print('Integration methods adapted')
md.solve()
U = md.variable('u')
print("iter %g"%(i+1))
_______________________________________________
Getfem-users mailing list
[email protected]
https://mail.gna.org/listinfo/getfem-users
