See for instance: Stokes equations with Mini elements — DOLFIN documentation
Where the MeshFunction sub_domains is used to apply boundary conditions.
For using meshfunctions in integration measures, see: Two ways to use MeshFunction to mark domain boundaries - #2 by dokken
I tried with a simpler mesh found on internet and the program test colab works. I think i didn’t correctly defined my mesh on gmsh. I am a student and i am new to gmsh. Could you just help me to fix my gmsh script ?
Jt = 0.2 ;
Ht = 0.2 ;
D = 1. ;
B = 3*D ;
L = 5*D ;
S = D ;
H = Ht*D ;
jmax = Sqrt( ((D/2)^2)+H^2) ;
j = Jt*jmax ;
xI = -(((j/2)/jmax)*D)/2 ;
yI = -Sqrt(((j/2)^2)-(xI^2)) ;
xJ = -xI ;
yJ = -yI ;
xF = xI-Sqrt(((jmax-j)^2)-(H-2*yJ)^2) ;
yF = yI-(H-2*yJ) ;
xM = xJ + Sqrt(((jmax-j)^2)-(H-2*yJ)^2) ;
yM = yJ + (H-2*yJ) ;
xL = xM-D ;
yL = yM ;
xG = xF+D ;
yG = yF ;
xE = -(L/2) ;
yE = yF ;
xH = L/2 ;
yH = yF ;
xK = xE ;
yK = yL ;
xN = xH ;
yN = yL ;
xC = xE ;
yC = yE-S ;
xD = xH ;
yD = yH-S ;
xO = xK ;
yO = yK+S ;
xP = xN ;
yP = yN+S ;
xA = xC ;
yA = yC-(B-S) ;
xB = xD ;
yB = yD-(B-S) ;
xQ = xO ;
yQ = yO+(B-S) ;
xR = xP ;
yR = yP+(B-S) ;
Point(1) = {xA, yA, 0, 1};
Point(2) = {xB, yB, 0, 1};
Point(3) = {xC, yC, 0, 1};
Point(4) = {xD, yD, 0, 1};
Point(5) = {xE, yE, 0, 1};
Point(6) = {xF, yF, 0, 1};
Point(7) = {xG, yG, 0, 1};
Point(8) = {xH, yH, 0, 1};
Point(9) = {xI, yI, 0, 1};
Point(10) = {xJ, yJ, 0, 1};
Point(11) = {xK, yK, 0, 1};
Point(12) = {xL, yL, 0, 1};
Point(13) = {xM, yM, 0, 1};
Point(14) = {xN, yN, 0, 1};
Point(15) = {xO, yO, 0, 1};
Point(16) = {xP, yP, 0, 1};
Point(17) = {xQ, yQ, 0, 1};
Point(18) = {xR, yR, 0, 1};
Line(1) = {1, 2};
Line(2) = {2, 4};
Line(3) = {8, 7};
Line(4) = {7, 10};
Line(5) = {10, 6};
Line(6) = {6, 5};
Line(7) = {3, 1};
Line(8) = {9, 12};
Line(9) = {12, 11};
Line(10) = {15, 17};
Line(11) = {17, 18};
Line(12) = {18, 16};
Line(13) = {14, 13};
Line(14) = {13, 9};
Line(15) = {15, 16};
Line(16) = {4, 3};
Line(17) = {11, 15};
Line(18) = {14, 16};
Line(19) = {5, 3};
Line(20) = {4, 8};
Line(21) = {12, 13};
Line(22) = {7, 6};
Line(55) = {10, 13};
Line(144) = {9, 6};
Line Loop(23) = {1, 2, 16, 7};
Line Loop(24) = {-16, 20, 3, 22, 6, 19};
Line Loop(25) = {-22, 4, 5};
Line Loop(26) = {21, 8, 14};
Line Loop(27) = {-15, -17, -9, 21, -13, 18};
Line Loop(28) = {15, -12, -11, -10};
Plane Surface(29) = {23};
Plane Surface(30) = {24};
Plane Surface(31) = {25};
Plane Surface(32) = {26};
Plane Surface(33) = {27};
Plane Surface(34) = {28};
Transfinite Line {1, 11} = 5 ;
Transfinite Line {-12, 10, 7, -2} = 2*5 Using Progression 1;
Transfinite Line {15, 16} = 500 ;
Transfinite Line {17, 18, 19, 20, 21, 22, 8, 14, 5, 4} = 100;
Transfinite Line {9, 13, 6, 3} = 200;
Physical Line(35) = {1} ;
Physical Line(36) = {2} ;
Physical Line(37) = {20} ;
Physical Line(38) = {3} ;
Physical Line(39) = {6} ;
Physical Line(40) = {19} ;
Physical Line(41) = {7} ;
Physical Line(42) = {10} ;
Physical Line(43) = {11} ;
Physical Line(44) = {12} ;
Physical Line(45) = {18} ;
Physical Line(46) = {13} ;
Physical Line(47) = {9} ;
Physical Line(48) = {17} ;
Physical Line(49) = {8} ;
Physical Line(50) = {4} ;
Physical Line(51) = {55} ;
Physical Line(52) = {144} ;
Physical Surface(53) = {29};
Physical Surface(54) = {30};
Physical Surface(56) = {31};
Physical Surface(57) = {32};
Physical Surface(58) = {33};
Physical Surface(59) = {34};
The geometry is this one, the bc are prescribed displacements and are in green
Jt and Ht are parameters that i want to control further in order to consider different meshes. They control the geometry of the triangles Jt=J/Jmax and Ht = D/H…
And D is the distance between F and G. It is the same between L and M
And this is how i would like to solve the problem, but i have to define “top” and “bottom” correctly, linked to the boundaries of the gmsh mesh (green parts)
# define elastic problem
dEu = derivative(total_en, u, v)
dEuu = derivative(dEu, u, u_)
# boundary conditions in displacement
bottom = CompiledSubDomain("near(x[1], %s, DOLFIN_EPS)"%(yF-B))
top = CompiledSubDomain("near(x[1], %s, DOLFIN_EPS)"%(yL+B))
g = Expression(("t", 0.), t=1., degree=0)
bottombound = project(-g, V_u)
topbound = project(g, V_u)
bc_u0 = DirichletBC(V_u, -g, bottom)
bc_u1 = DirichletBC(V_u, g, top)
bc_u = [bc_u0, bc_u1]
# solver for the elastic problem
problem_u = NonlinearVariationalProblem(dEu, u, bc_u, dEuu)
solver_u = NonlinearVariationalSolver(problem_u)
solver_u.parameters.update({"newton_solver":{ "linear_solver" : "mumps"}})
# define damage problem
dEa = derivative(total_en, alpha, beta)
dEaa = derivative(dEa, alpha, dalpha)
# boundary conditions in damage
upperbound = project(Constant(1.), V_alpha)
lowerbound = project(Constant(0.), V_alpha)
bc_a1 = DirichletBC(V_alpha, Constant(0.), bottom)
bc_a2 = DirichletBC(V_alpha, Constant(0.), top)
bc_a = [bc_a1, bc_a2]
# solver for the damage problem
problem_alpha = NonlinearVariationalProblem(dEa, alpha, bc_a, dEaa)
problem_alpha.set_bounds(lowerbound, upperbound)
solver_alpha = NonlinearVariationalSolver(problem_alpha)
solver_alpha.parameters.update(params_alpha)