I'm trying to generate a 3D mesh using the python API of gmsh, and then use this mesh on a FEM software (Code_Aster). My purpose here is to create an articifical layer (a joint) between 2 volumes (soil & foundation) to implement cohesive elements in my model. My problem is that it seems like in my FEM software, I need the elements from my joint volume to be of PENTA6 type (prism). And in python gmsh, the default elements created are tetrahedral. Thus my question is how could I change the element type for the joint volume so that they become prisms and work on my FEM model?
I tried performing extrusions and using the "recombine" keyword to create my joint, and I managed to create prism elements but my FEM software seems to have problems understanding these extrusions.
ext_jointz1 = gmsh.model.geo.extrude(dimTags_jointz1, 0, 0, e_joint, [1], recombine=True)
My code as of right now:
import gmsh
import math
import sys
import numpy as np
# Before using any functions in the Python API, Gmsh must be initialized:
gmsh.initialize()
###### Geometrical Constants #######
# SI units
Supercell = 10 # "Sol" dimensions
Ox = Supercell/2 # Foundation center x-coordinate
Oy = Supercell/2 # Foundation center y-coordinate
# Foundation dimensions given by RTE
a = 1.7
b = 2.1
c = 2.1
d = 2.9
e = 2.9
f = 0.4
h = 0.3 # the foundation rises out of the ground by 30cm
delta = 0.25 # for measuring points
HTube = 4 # height of the tube
e_joint = 0.01 # thickness of the joint layer
###### Mesh Constants #######
lc_joint = 0.1
lc_prec = 0.1
lc = 0.5
lc_i = 1
lc_out = 2
###### Geometry Definition #######
gmsh.clear()
## Add points
#-------------------------------------#
#------- domain center --------------#
#-------------------------------------#
# gmsh.model.geo.addPoint(0, 0, 0, lc, 0)
# Soil points
p1 = gmsh.model.geo.addPoint(0, 0, 0, lc_out, 1);
p2 = gmsh.model.geo.addPoint(Supercell, 0, 0, lc_out, 2);
p3 = gmsh.model.geo.addPoint(Supercell, Supercell, 0, lc_out, 3);
p4 = gmsh.model.geo.addPoint(0, Supercell, 0, lc_out, 4);
p5 = gmsh.model.geo.addPoint(0, 0, Supercell, lc_out, 5);
p6 = gmsh.model.geo.addPoint(Supercell, 0, Supercell, lc_out, 6);
p7 = gmsh.model.geo.addPoint(Supercell, Supercell, Supercell, lc_out, 7);
p8 = gmsh.model.geo.addPoint(0, Supercell, Supercell, lc_out, 8);
# Foundation points
# Socle (partie basse de la fondation)
# Points du bas du socle
p9 = gmsh.model.geo.addPoint(Ox+e/2, Oy-d/2, Supercell-a, lc_prec, 9);
p10 = gmsh.model.geo.addPoint(Ox+e/2, Oy+d/2, Supercell-a, lc_prec, 10);
p11 = gmsh.model.geo.addPoint(Ox-e/2, Oy+d/2, Supercell-a, lc_prec, 11);
p12 = gmsh.model.geo.addPoint(Ox-e/2, Oy-d/2, Supercell-a, lc_prec, 12);
# Points du haut du socle
p13 = gmsh.model.geo.addPoint(Ox+e/2, Oy-d/2, Supercell-a+f, lc_prec, 13);
p14 = gmsh.model.geo.addPoint(Ox+e/2, Oy+d/2, Supercell-a+f, lc_prec, 14);
p15 = gmsh.model.geo.addPoint(Ox-e/2, Oy+d/2, Supercell-a+f, lc_prec, 15);
p16 = gmsh.model.geo.addPoint(Ox-e/2, Oy-d/2, Supercell-a+f, lc_prec, 16);
# Base supérieure (partie haute de la fondation)
p17 = gmsh.model.geo.addPoint(Ox+c/2, Oy-b/2, Supercell-a+f, lc_prec, 17);
p18 = gmsh.model.geo.addPoint(Ox+c/2, Oy+b/2, Supercell-a+f, lc_prec, 18);
p19 = gmsh.model.geo.addPoint(Ox-c/2, Oy+b/2, Supercell-a+f, lc_prec, 19);
p20 = gmsh.model.geo.addPoint(Ox-c/2, Oy-b/2, Supercell-a+f, lc_prec, 20);
p21 = gmsh.model.geo.addPoint(Ox+c/2, Oy-b/2, Supercell, lc_prec, 21);
p22 = gmsh.model.geo.addPoint(Ox+c/2, Oy+b/2, Supercell, lc_prec, 22);
p23 = gmsh.model.geo.addPoint(Ox-c/2, Oy+b/2, Supercell, lc_prec, 23);
p24 = gmsh.model.geo.addPoint(Ox-c/2, Oy-b/2, Supercell, lc_prec, 24);
p25 = gmsh.model.geo.addPoint(Ox+c/2, Oy-b/2, Supercell+h, lc, 25);
p26 = gmsh.model.geo.addPoint(Ox+c/2, Oy+b/2, Supercell+h, lc, 26);
p27 = gmsh.model.geo.addPoint(Ox-c/2, Oy+b/2, Supercell+h, lc, 27);
p28 = gmsh.model.geo.addPoint(Ox-c/2, Oy-b/2, Supercell+h, lc, 28);
# Measuring points
p29 = gmsh.model.geo.addPoint(Ox+1, Oy+1, Supercell+h, lc_prec, 29);
p30 = gmsh.model.geo.addPoint(Ox-1, Oy+1, Supercell+h, lc_prec, 30);
p31 = gmsh.model.geo.addPoint(Ox-1, Oy-1, Supercell+h, lc_prec, 31);
p32 = gmsh.model.geo.addPoint(Ox+1, Oy-1, Supercell+h, lc_prec, 32);
p33 = gmsh.model.geo.addPoint(Ox+1-delta, Oy+1-delta, Supercell+h, lc_prec, 33);
p34 = gmsh.model.geo.addPoint(Ox-1+delta, Oy+1-delta, Supercell+h, lc_prec, 34);
p35 = gmsh.model.geo.addPoint(Ox-1+delta, Oy-1+delta, Supercell+h, lc_prec, 35);
p36 = gmsh.model.geo.addPoint(Ox+1-delta, Oy-1+delta, Supercell+h, lc_prec, 36);
# Center point
p37 = gmsh.model.geo.addPoint(Ox, Oy, Supercell+h, lc, 37);
# Tube point (remote force point of application)
p38 = gmsh.model.geo.addPoint(Ox, Oy, Supercell+h+HTube, lc, 38);
# to be checked: we double the points of the HOLE, defining contact between the two volumes/physicals?
# Joint points
# Socle (partie basse de la fondation)
# Points du bas du socle
p39 = gmsh.model.geo.addPoint(Ox+e/2+e_joint, Oy-d/2-e_joint, Supercell-a-e_joint, lc_joint, 39);
p40 = gmsh.model.geo.addPoint(Ox+e/2+e_joint, Oy+d/2+e_joint, Supercell-a-e_joint, lc_joint, 40);
p41 = gmsh.model.geo.addPoint(Ox-e/2-e_joint, Oy+d/2+e_joint, Supercell-a-e_joint, lc_joint, 41);
p42 = gmsh.model.geo.addPoint(Ox-e/2-e_joint, Oy-d/2-e_joint, Supercell-a-e_joint, lc_joint, 42);
# Points du haut du socle
p43 = gmsh.model.geo.addPoint(Ox+e/2+e_joint, Oy-d/2-e_joint, Supercell-a+f+e_joint, lc_joint, 43);
p44 = gmsh.model.geo.addPoint(Ox+e/2+e_joint, Oy+d/2+e_joint, Supercell-a+f+e_joint, lc_joint, 44);
p45 = gmsh.model.geo.addPoint(Ox-e/2-e_joint, Oy+d/2+e_joint, Supercell-a+f+e_joint, lc_joint, 45);
p46 = gmsh.model.geo.addPoint(Ox-e/2-e_joint, Oy-d/2-e_joint, Supercell-a+f+e_joint, lc_joint, 46);
# Base supérieure (partie haute de la fondation)
p47 = gmsh.model.geo.addPoint(Ox+c/2+e_joint, Oy-b/2-e_joint, Supercell-a+f+e_joint, lc_joint, 47);
p48 = gmsh.model.geo.addPoint(Ox+c/2+e_joint, Oy+b/2+e_joint, Supercell-a+f+e_joint, lc_joint, 48);
p49 = gmsh.model.geo.addPoint(Ox-c/2-e_joint, Oy+b/2+e_joint, Supercell-a+f+e_joint, lc_joint, 49);
p50 = gmsh.model.geo.addPoint(Ox-c/2-e_joint, Oy-b/2-e_joint, Supercell-a+f+e_joint, lc_joint, 50);
p51 = gmsh.model.geo.addPoint(Ox+c/2+e_joint, Oy-b/2-e_joint, Supercell, lc_joint, 51);
p52 = gmsh.model.geo.addPoint(Ox+c/2+e_joint, Oy+b/2+e_joint, Supercell, lc_joint, 52);
p53 = gmsh.model.geo.addPoint(Ox-c/2-e_joint, Oy+b/2+e_joint, Supercell, lc_joint, 53);
p54 = gmsh.model.geo.addPoint(Ox-c/2-e_joint, Oy-b/2-e_joint, Supercell, lc_joint, 54);
## Add lines
# Soil lines
l1 = gmsh.model.geo.addLine(1, 2, 1)
l2 = gmsh.model.geo.addLine(2, 3, 2)
l3 = gmsh.model.geo.addLine(3, 4, 3)
l4 = gmsh.model.geo.addLine(4, 1, 4)
l5 = gmsh.model.geo.addLine(1, 5, 5)
l6 = gmsh.model.geo.addLine(2, 6, 6)
l7 = gmsh.model.geo.addLine(3, 7, 7)
l8 = gmsh.model.geo.addLine(4, 8, 8)
l9 = gmsh.model.geo.addLine(5, 6, 9)
l10 = gmsh.model.geo.addLine(6, 7, 10)
l11 = gmsh.model.geo.addLine(7, 8, 11)
l12 = gmsh.model.geo.addLine(8, 5, 12)
# Foundation lines
l13 = gmsh.model.geo.addLine(9, 10, 13)
l14 = gmsh.model.geo.addLine(10, 11, 14)
l15 = gmsh.model.geo.addLine(11, 12, 15)
l16 = gmsh.model.geo.addLine(12, 9, 16)
l17 = gmsh.model.geo.addLine(9, 13, 17)
l18 = gmsh.model.geo.addLine(10, 14, 18)
l19 = gmsh.model.geo.addLine(11, 15, 19)
l20 = gmsh.model.geo.addLine(12, 16, 20)
l21 = gmsh.model.geo.addLine(13, 14, 21)
l22 = gmsh.model.geo.addLine(14, 15, 22)
l23 = gmsh.model.geo.addLine(15, 16, 23)
l24 = gmsh.model.geo.addLine(16, 13, 24)
l25 = gmsh.model.geo.addLine(17, 18, 25)
l26 = gmsh.model.geo.addLine(18, 19, 26)
l27 = gmsh.model.geo.addLine(19, 20, 27)
l28 = gmsh.model.geo.addLine(20, 17, 28)
l29 = gmsh.model.geo.addLine(17, 21, 29)
l30 = gmsh.model.geo.addLine(18, 22, 30)
l31 = gmsh.model.geo.addLine(19, 23, 31)
l32 = gmsh.model.geo.addLine(20, 24, 32)
l33 = gmsh.model.geo.addLine(21, 22, 33)
l34 = gmsh.model.geo.addLine(22, 23, 34)
l35 = gmsh.model.geo.addLine(23, 24, 35)
l36 = gmsh.model.geo.addLine(24, 21, 36)
l37 = gmsh.model.geo.addLine(21, 25, 37)
l38 = gmsh.model.geo.addLine(22, 26, 38)
l39 = gmsh.model.geo.addLine(23, 27, 39)
l40 = gmsh.model.geo.addLine(24, 28, 40)
l41 = gmsh.model.geo.addLine(25, 26, 41)
l42 = gmsh.model.geo.addLine(26, 27, 42)
l43 = gmsh.model.geo.addLine(27, 28, 43)
l44 = gmsh.model.geo.addLine(28, 25, 44)
# Joint lines
l45 = gmsh.model.geo.addLine(39, 40, 45)
l46 = gmsh.model.geo.addLine(40, 41, 46)
l47 = gmsh.model.geo.addLine(41, 42, 47)
l48 = gmsh.model.geo.addLine(42, 39, 48)
l49 = gmsh.model.geo.addLine(39, 43, 49)
l50 = gmsh.model.geo.addLine(40, 44, 50)
l51 = gmsh.model.geo.addLine(41, 45, 51)
l52 = gmsh.model.geo.addLine(42, 46, 52)
l53 = gmsh.model.geo.addLine(43, 44, 53)
l54 = gmsh.model.geo.addLine(44, 45, 54)
l55 = gmsh.model.geo.addLine(45, 46, 55)
l56 = gmsh.model.geo.addLine(46, 43, 56)
l57 = gmsh.model.geo.addLine(47, 48, 57)
l58 = gmsh.model.geo.addLine(48, 49, 58)
l59 = gmsh.model.geo.addLine(49, 50, 59)
l60 = gmsh.model.geo.addLine(50, 47, 60)
l61 = gmsh.model.geo.addLine(47, 51, 61)
l62 = gmsh.model.geo.addLine(48, 52, 62)
l63 = gmsh.model.geo.addLine(49, 53, 63)
l64 = gmsh.model.geo.addLine(50, 54, 64)
l65 = gmsh.model.geo.addLine(51, 52, 65)
l66 = gmsh.model.geo.addLine(52, 53, 66)
l67 = gmsh.model.geo.addLine(53, 54, 67)
l68 = gmsh.model.geo.addLine(54, 51, 68)
## Add Curve Loops and Plane Surfaces
# Soil #
# Soil bottom plane
cl1 = gmsh.model.geo.addCurveLoop([1, 2, 3, 4])
s1 = gmsh.model.geo.addPlaneSurface([cl1])
# Soil side planes
cl2 = gmsh.model.geo.addCurveLoop([1, 6, -9, -5])
s2 = gmsh.model.geo.addPlaneSurface([cl2])
cl3 = gmsh.model.geo.addCurveLoop([2, 7, -10, -6])
s3 = gmsh.model.geo.addPlaneSurface([cl3])
cl4 = gmsh.model.geo.addCurveLoop([3, 8, -11, -7])
s4 = gmsh.model.geo.addPlaneSurface([cl4])
cl5 = gmsh.model.geo.addCurveLoop([4, 5, -12, -8])
s5 = gmsh.model.geo.addPlaneSurface([cl5])
# Soil upper plane
cl600 = gmsh.model.geo.addCurveLoop([9, 10, 11, 12])
cl60 = gmsh.model.geo.addCurveLoop([65, 66, 67, 68])
s6 = gmsh.model.geo.addPlaneSurface([cl600, cl60])
# Foundation #
# Foundation bottom side planes
cl7 = gmsh.model.geo.addCurveLoop([13, 18, -21, -17])
s7 = gmsh.model.geo.addPlaneSurface([cl7])
cl8 = gmsh.model.geo.addCurveLoop([14, 19, -22, -18])
s8 = gmsh.model.geo.addPlaneSurface([cl8])
cl9 = gmsh.model.geo.addCurveLoop([15, 20, -23, -19])
s9 = gmsh.model.geo.addPlaneSurface([cl9])
cl10 = gmsh.model.geo.addCurveLoop([16, 17, -24, -20])
s10 = gmsh.model.geo.addPlaneSurface([cl10])
# Foundation bottom plane
cl11 = gmsh.model.geo.addCurveLoop([13, 14, 15, 16])
s11 = gmsh.model.geo.addPlaneSurface([cl11])
# Foundation bottom upper plane
cl1200 = gmsh.model.geo.addCurveLoop([21, 22, 23, 24])
cl120 = gmsh.model.geo.addCurveLoop([25, 26, 27, 28])
s12 = gmsh.model.geo.addPlaneSurface([cl1200, cl120])
# Foundation upper bottom side planes
cl13 = gmsh.model.geo.addCurveLoop([25, 30, -33, -29])
s13 = gmsh.model.geo.addPlaneSurface([cl13])
cl14 = gmsh.model.geo.addCurveLoop([26, 31, -34, -30])
s14 = gmsh.model.geo.addPlaneSurface([cl14])
cl15 = gmsh.model.geo.addCurveLoop([27, 32, -35, -31])
s15 = gmsh.model.geo.addPlaneSurface([cl15])
cl16 = gmsh.model.geo.addCurveLoop([28, 29, -36, -32])
s16 = gmsh.model.geo.addPlaneSurface([cl16])
# Foundation upper upper side planes
cl17 = gmsh.model.geo.addCurveLoop([33, 38, -41, -37])
s17 = gmsh.model.geo.addPlaneSurface([cl17])
cl18 = gmsh.model.geo.addCurveLoop([34, 39, -42, -38])
s18 = gmsh.model.geo.addPlaneSurface([cl18])
cl19 = gmsh.model.geo.addCurveLoop([35, 40, -43, -39])
s19 = gmsh.model.geo.addPlaneSurface([cl19])
cl20 = gmsh.model.geo.addCurveLoop([36, 37, -44, -40])
s20 = gmsh.model.geo.addPlaneSurface([cl20])
# Foundation upper plane
cl21 = gmsh.model.geo.addCurveLoop([41, 42, 43, 44])
s21 = gmsh.model.geo.addPlaneSurface([cl21])
gmsh.model.geo.synchronize()
gmsh.model.mesh.embed(0, [29, 30, 31, 32, 33, 34, 35, 36], 2, 21) # add measuring points to upper surface of the foundation
# Joint #
# Joint bottom side planes
cl22 = gmsh.model.geo.addCurveLoop([45, 50, -53, -49])
s22 = gmsh.model.geo.addPlaneSurface([cl22])
cl23 = gmsh.model.geo.addCurveLoop([46, 51, -54, -50])
s23 = gmsh.model.geo.addPlaneSurface([cl23])
cl24 = gmsh.model.geo.addCurveLoop([47, 52, -55, -51])
s24 = gmsh.model.geo.addPlaneSurface([cl24])
cl25 = gmsh.model.geo.addCurveLoop([48, 49, -56, -52])
s25 = gmsh.model.geo.addPlaneSurface([cl25])
# Joint bottom plane
cl26 = gmsh.model.geo.addCurveLoop([45, 46, 47, 48])
s26 = gmsh.model.geo.addPlaneSurface([cl26])
# Joint bottom upper plane
cl2700 = gmsh.model.geo.addCurveLoop([53, 54, 55, 56])
cl270 = gmsh.model.geo.addCurveLoop([57, 58, 59, 60])
s27 = gmsh.model.geo.addPlaneSurface([cl2700, cl270])
# Joint upper bottom side planes
cl28 = gmsh.model.geo.addCurveLoop([57, 62, -65, -61])
s28 = gmsh.model.geo.addPlaneSurface([cl28])
cl29 = gmsh.model.geo.addCurveLoop([58, 63, -66, -62])
s29 = gmsh.model.geo.addPlaneSurface([cl29])
cl30 = gmsh.model.geo.addCurveLoop([59, 64, -67, -63])
s30 = gmsh.model.geo.addPlaneSurface([cl30])
cl31 = gmsh.model.geo.addCurveLoop([60, 61, -68, -64])
s31 = gmsh.model.geo.addPlaneSurface([cl31])
# Joint upper plane
cl3200 = gmsh.model.geo.addCurveLoop([65, 66, 67, 68])
cl320 = gmsh.model.geo.addCurveLoop([33, 34, 35, 36])
s32 = gmsh.model.geo.addPlaneSurface([cl3200, cl320])
gmsh.model.geo.synchronize()
## Add Volume
# Soil volume
sl1 = gmsh.model.geo.addSurfaceLoop([1, 2, 3, 4, 5, 6, 28, 29, 30, 31, 27, 22, 23, 24, 25, 26])
v1 = gmsh.model.geo.addVolume([sl1])
# Foundation volume
sl2 = gmsh.model.geo.addSurfaceLoop([11, 7, 8, 9, 10, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21])
v2 = gmsh.model.geo.addVolume([sl2])
# Joint volume
# gmsh.option.setNumber("Mesh.Algorithm", 2);
sl3 = gmsh.model.geo.addSurfaceLoop([26, 22, 23, 24, 25, 27, 28, 29, 30, 31, 32, 13, 14, 15, 16, 12, 7, 8, 9, 10, 11])
v3 = gmsh.model.geo.addVolume([sl3])
## Add Physical Groups
# Soil
pv1 = gmsh.model.addPhysicalGroup(3, [1])
gmsh.model.setPhysicalName(3, pv1, "SOIL")
# Foundation
pv2 = gmsh.model.addPhysicalGroup(3, [2])
gmsh.model.setPhysicalName(3, pv2, "MASSIF")
# Joint
pv3 = gmsh.model.addPhysicalGroup(3, [3])
gmsh.model.setPhysicalName(3, pv3, "JOINT")
# Soil bottom
ps1 = gmsh.model.addPhysicalGroup(2, [1])
gmsh.model.setPhysicalName(2, ps1, "BAS")
# Upper Soil
ps2 = gmsh.model.addPhysicalGroup(2, [6])
gmsh.model.setPhysicalName(2, ps2, "UPSOL")
# Upper Foundation
ps3 = gmsh.model.addPhysicalGroup(2, [21])
gmsh.model.setPhysicalName(2, ps3, "UPMASSIF")
# Upper Joint
ps4 = gmsh.model.addPhysicalGroup(2, [32])
gmsh.model.setPhysicalName(2, ps4, "UPJOINT")
# Soil EXT
ps5 = gmsh.model.addPhysicalGroup(2, [2, 3, 4, 5])
gmsh.model.setPhysicalName(2, ps5, "EXT")
# Contact Soil/Joint
ps6 = gmsh.model.addPhysicalGroup(2, [22, 23, 24, 25, 26, 27, 28, 29, 30, 31])
gmsh.model.setPhysicalName(2, ps6, "CONTACT_SJ")
# Contact Joint/Massif
ps7 = gmsh.model.addPhysicalGroup(2, [7, 8, 9, 10, 11, 12, 13, 14, 15, 16])
gmsh.model.setPhysicalName(2, ps7, "CONTACT_JM")
# ?
pp1 = gmsh.model.addPhysicalGroup(0, [9])
gmsh.model.setPhysicalName(0, pp1, "APPLI")
# Measuring points
pp2 = gmsh.model.addPhysicalGroup(0, [29, 30, 31, 32])
gmsh.model.setPhysicalName(0, pp2, "A")
pp3 = gmsh.model.addPhysicalGroup(0, [33, 34, 35, 36])
gmsh.model.setPhysicalName(0, pp3, "Abis")
# Center point
pp4 = gmsh.model.addPhysicalGroup(0, [37])
gmsh.model.setPhysicalName(0, pp4, "Center")
# Remote force application point
pp5 = gmsh.model.addPhysicalGroup(0, [38])
gmsh.model.setPhysicalName(0, pp5, "Tube")
# gmsh.model.geo.removeAllDuplicates()
# gmsh.model.mesh.setRecombine(3, 3)
######## Mesh Generation #############
# gmsh.option.setNumber("Mesh.RecombineAll", 1)
gmsh.model.geo.synchronize()
gmsh.model.mesh.generate()
gmsh.option.setNumber("Mesh.MedFileMinorVersion", 0)
gmsh.write("FE_model_mesh51.med")
# gmsh.model.mesh.removeDuplicateNodes()
# gmsh.write("../mesh_per_angle_refined/refined_angle="+list4[list_alphap.index(alphap)]+".msh")
############################
gmsh.model.geo.synchronize()
#gmsh.fltk.run()
gmsh.finalize()