Skip to content
Snippets Groups Projects
Commit 6e1eea89 authored by Mohib's avatar Mohib
Browse files

[Minor] - HouseKeeping

parent 7040c5a2
No related branches found
No related tags found
1 merge request!433[GenericCrackPhaseField] - Parity with new pf interface
Showing
with 0 additions and 41469 deletions
//Fiber reinforced composite cross-section
Geometry.AutoCoherence= 1.0e-9;
//defination of unit
unit = 1.0e-3;
//characteristic mesh size
h_x=0.05*unit;
h_y=0.06*unit;//0.1
h_z = 0.02*unit;
L = 0.7*unit;
nx= 1;
ny= 1;
cl=0.1*unit;
//point
Point(1) = {0.0 , 0.0 , 0.0, cl};
Point(2) = {0.0 , ny*h_y , 0.0, cl};
Line(1) = {1,2};
Extrude {nx*h_x , 0.0 , 0.0} {Line{1}; Layers{nx}; Recombine;}
Extrude {L*2 , 0.0 , 0.0} {Line{2};Layers{60}; Recombine;}
Extrude {nx*h_x , 0.0 , 0.0} {Line{6}; Layers{nx}; Recombine;}
Physical Line(100) = {1};
Physical Line(101) = {10};
Physical Line(110) = {3,7,11};
Physical Line(111) = {4,8,12};
Physical Surface(120) = {5,9,13};
Physical Point(1011) = {1};
Physical Point(1012) = {2};
Physical Surface(1041) = {5,13};
Physical Point(1015) = {15};
Physical Point(1016) = {16};
Coherence;
$MeshFormat
3 0 8
$EndMeshFormat
$Entities
8 10 3 0
1 1 1011
2 1 1012
3 0
4 0
9 0
10 0
15 1 1015
16 1 1016
1 2 1 2 1 100
2 2 3 4 0
3 2 1 3 1 110
4 2 2 4 1 111
6 2 9 10 0
7 2 3 9 1 110
8 2 4 10 1 111
10 2 15 16 1 101
11 2 9 15 1 110
12 2 10 16 1 111
5 4 1 4 -2 -3 2 120 1041
9 4 2 8 -6 -7 1 120
13 4 6 12 -10 -11 2 120 1041
$EndEntities
$Nodes
126
1 0 0 0 0
2 0 6e-05 0 0
3 5e-05 0 0 0
4 5e-05 6e-05 0 0
5 0.00145 0 0 0
6 0.00145 6e-05 0 0
7 0.0015 0 0 0
8 0.0015 6e-05 0 0
9 7.333333333333333e-05 0 0 0
10 9.666666666666667e-05 0 0 0
11 0.00012 0 0 0
12 0.0001433333333333333 0 0 0
13 0.0001666666666666667 0 0 0
14 0.00019 0 0 0
15 0.0002133333333333333 0 0 0
16 0.0002366666666666667 0 0 0
17 0.00026 0 0 0
18 0.0002833333333333334 0 0 0
19 0.0003066666666666667 0 0 0
20 0.0003300000000000001 0 0 0
21 0.0003533333333333334 0 0 0
22 0.0003766666666666667 0 0 0
23 0.0004 0 0 0
24 0.0004233333333333333 0 0 0
25 0.0004466666666666667 0 0 0
26 0.00047 0 0 0
27 0.0004933333333333333 0 0 0
28 0.0005166666666666667 0 0 0
29 0.00054 0 0 0
30 0.0005633333333333333 0 0 0
31 0.0005866666666666668 0 0 0
32 0.0006100000000000001 0 0 0
33 0.0006333333333333334 0 0 0
34 0.0006566666666666667 0 0 0
35 0.00068 0 0 0
36 0.0007033333333333334 0 0 0
37 0.0007266666666666667 0 0 0
38 0.00075 0 0 0
39 0.0007733333333333334 0 0 0
40 0.0007966666666666667 0 0 0
41 0.0008200000000000001 0 0 0
42 0.0008433333333333333 0 0 0
43 0.0008666666666666667 0 0 0
44 0.0008899999999999999 0 0 0
45 0.0009133333333333334 0 0 0
46 0.0009366666666666666 0 0 0
47 0.00096 0 0 0
48 0.0009833333333333332 0 0 0
49 0.001006666666666667 0 0 0
50 0.00103 0 0 0
51 0.001053333333333333 0 0 0
52 0.001076666666666667 0 0 0
53 0.0011 0 0 0
54 0.001123333333333333 0 0 0
55 0.001146666666666666 0 0 0
56 0.00117 0 0 0
57 0.001193333333333333 0 0 0
58 0.001216666666666667 0 0 0
59 0.00124 0 0 0
60 0.001263333333333333 0 0 0
61 0.001286666666666666 0 0 0
62 0.00131 0 0 0
63 0.001333333333333333 0 0 0
64 0.001356666666666667 0 0 0
65 0.00138 0 0 0
66 0.001403333333333333 0 0 0
67 0.001426666666666667 0 0 0
68 7.333333333333333e-05 6e-05 0 0
69 9.666666666666667e-05 6e-05 0 0
70 0.00012 6e-05 0 0
71 0.0001433333333333333 6e-05 0 0
72 0.0001666666666666667 6e-05 0 0
73 0.00019 6e-05 0 0
74 0.0002133333333333333 6e-05 0 0
75 0.0002366666666666667 6e-05 0 0
76 0.00026 6e-05 0 0
77 0.0002833333333333334 6e-05 0 0
78 0.0003066666666666667 6e-05 0 0
79 0.0003300000000000001 6e-05 0 0
80 0.0003533333333333334 6e-05 0 0
81 0.0003766666666666667 6e-05 0 0
82 0.0004 6e-05 0 0
83 0.0004233333333333333 6e-05 0 0
84 0.0004466666666666667 6e-05 0 0
85 0.00047 6e-05 0 0
86 0.0004933333333333333 6e-05 0 0
87 0.0005166666666666667 6e-05 0 0
88 0.00054 6e-05 0 0
89 0.0005633333333333333 6e-05 0 0
90 0.0005866666666666668 6e-05 0 0
91 0.0006100000000000001 6e-05 0 0
92 0.0006333333333333334 6e-05 0 0
93 0.0006566666666666667 6e-05 0 0
94 0.00068 6e-05 0 0
95 0.0007033333333333334 6e-05 0 0
96 0.0007266666666666667 6e-05 0 0
97 0.00075 6e-05 0 0
98 0.0007733333333333334 6e-05 0 0
99 0.0007966666666666667 6e-05 0 0
100 0.0008200000000000001 6e-05 0 0
101 0.0008433333333333333 6e-05 0 0
102 0.0008666666666666667 6e-05 0 0
103 0.0008899999999999999 6e-05 0 0
104 0.0009133333333333334 6e-05 0 0
105 0.0009366666666666666 6e-05 0 0
106 0.00096 6e-05 0 0
107 0.0009833333333333332 6e-05 0 0
108 0.001006666666666667 6e-05 0 0
109 0.00103 6e-05 0 0
110 0.001053333333333333 6e-05 0 0
111 0.001076666666666667 6e-05 0 0
112 0.0011 6e-05 0 0
113 0.001123333333333333 6e-05 0 0
114 0.001146666666666666 6e-05 0 0
115 0.00117 6e-05 0 0
116 0.001193333333333333 6e-05 0 0
117 0.001216666666666667 6e-05 0 0
118 0.00124 6e-05 0 0
119 0.001263333333333333 6e-05 0 0
120 0.001286666666666666 6e-05 0 0
121 0.00131 6e-05 0 0
122 0.001333333333333333 6e-05 0 0
123 0.001356666666666667 6e-05 0 0
124 0.00138 6e-05 0 0
125 0.001403333333333333 6e-05 0 0
126 0.001426666666666667 6e-05 0 0
$EndNodes
$Elements
192
131 3 5 4 1 2 4 3
132 3 9 4 3 4 68 9
133 3 9 4 9 68 69 10
134 3 9 4 10 69 70 11
135 3 9 4 11 70 71 12
136 3 9 4 12 71 72 13
137 3 9 4 13 72 73 14
138 3 9 4 14 73 74 15
139 3 9 4 15 74 75 16
140 3 9 4 16 75 76 17
141 3 9 4 17 76 77 18
142 3 9 4 18 77 78 19
143 3 9 4 19 78 79 20
144 3 9 4 20 79 80 21
145 3 9 4 21 80 81 22
146 3 9 4 22 81 82 23
147 3 9 4 23 82 83 24
148 3 9 4 24 83 84 25
149 3 9 4 25 84 85 26
150 3 9 4 26 85 86 27
151 3 9 4 27 86 87 28
152 3 9 4 28 87 88 29
153 3 9 4 29 88 89 30
154 3 9 4 30 89 90 31
155 3 9 4 31 90 91 32
156 3 9 4 32 91 92 33
157 3 9 4 33 92 93 34
158 3 9 4 34 93 94 35
159 3 9 4 35 94 95 36
160 3 9 4 36 95 96 37
161 3 9 4 37 96 97 38
162 3 9 4 38 97 98 39
163 3 9 4 39 98 99 40
164 3 9 4 40 99 100 41
165 3 9 4 41 100 101 42
166 3 9 4 42 101 102 43
167 3 9 4 43 102 103 44
168 3 9 4 44 103 104 45
169 3 9 4 45 104 105 46
170 3 9 4 46 105 106 47
171 3 9 4 47 106 107 48
172 3 9 4 48 107 108 49
173 3 9 4 49 108 109 50
174 3 9 4 50 109 110 51
175 3 9 4 51 110 111 52
176 3 9 4 52 111 112 53
177 3 9 4 53 112 113 54
178 3 9 4 54 113 114 55
179 3 9 4 55 114 115 56
180 3 9 4 56 115 116 57
181 3 9 4 57 116 117 58
182 3 9 4 58 117 118 59
183 3 9 4 59 118 119 60
184 3 9 4 60 119 120 61
185 3 9 4 61 120 121 62
186 3 9 4 62 121 122 63
187 3 9 4 63 122 123 64
188 3 9 4 64 123 124 65
189 3 9 4 65 124 125 66
190 3 9 4 66 125 126 67
191 3 9 4 67 126 6 5
192 3 13 4 5 6 8 7
5 1 1 2 1 2
6 1 3 2 1 3
7 1 4 2 2 4
8 1 7 2 3 9
9 1 7 2 9 10
10 1 7 2 10 11
11 1 7 2 11 12
12 1 7 2 12 13
13 1 7 2 13 14
14 1 7 2 14 15
15 1 7 2 15 16
16 1 7 2 16 17
17 1 7 2 17 18
18 1 7 2 18 19
19 1 7 2 19 20
20 1 7 2 20 21
21 1 7 2 21 22
22 1 7 2 22 23
23 1 7 2 23 24
24 1 7 2 24 25
25 1 7 2 25 26
26 1 7 2 26 27
27 1 7 2 27 28
28 1 7 2 28 29
29 1 7 2 29 30
30 1 7 2 30 31
31 1 7 2 31 32
32 1 7 2 32 33
33 1 7 2 33 34
34 1 7 2 34 35
35 1 7 2 35 36
36 1 7 2 36 37
37 1 7 2 37 38
38 1 7 2 38 39
39 1 7 2 39 40
40 1 7 2 40 41
41 1 7 2 41 42
42 1 7 2 42 43
43 1 7 2 43 44
44 1 7 2 44 45
45 1 7 2 45 46
46 1 7 2 46 47
47 1 7 2 47 48
48 1 7 2 48 49
49 1 7 2 49 50
50 1 7 2 50 51
51 1 7 2 51 52
52 1 7 2 52 53
53 1 7 2 53 54
54 1 7 2 54 55
55 1 7 2 55 56
56 1 7 2 56 57
57 1 7 2 57 58
58 1 7 2 58 59
59 1 7 2 59 60
60 1 7 2 60 61
61 1 7 2 61 62
62 1 7 2 62 63
63 1 7 2 63 64
64 1 7 2 64 65
65 1 7 2 65 66
66 1 7 2 66 67
67 1 7 2 67 5
68 1 8 2 4 68
69 1 8 2 68 69
70 1 8 2 69 70
71 1 8 2 70 71
72 1 8 2 71 72
73 1 8 2 72 73
74 1 8 2 73 74
75 1 8 2 74 75
76 1 8 2 75 76
77 1 8 2 76 77
78 1 8 2 77 78
79 1 8 2 78 79
80 1 8 2 79 80
81 1 8 2 80 81
82 1 8 2 81 82
83 1 8 2 82 83
84 1 8 2 83 84
85 1 8 2 84 85
86 1 8 2 85 86
87 1 8 2 86 87
88 1 8 2 87 88
89 1 8 2 88 89
90 1 8 2 89 90
91 1 8 2 90 91
92 1 8 2 91 92
93 1 8 2 92 93
94 1 8 2 93 94
95 1 8 2 94 95
96 1 8 2 95 96
97 1 8 2 96 97
98 1 8 2 97 98
99 1 8 2 98 99
100 1 8 2 99 100
101 1 8 2 100 101
102 1 8 2 101 102
103 1 8 2 102 103
104 1 8 2 103 104
105 1 8 2 104 105
106 1 8 2 105 106
107 1 8 2 106 107
108 1 8 2 107 108
109 1 8 2 108 109
110 1 8 2 109 110
111 1 8 2 110 111
112 1 8 2 111 112
113 1 8 2 112 113
114 1 8 2 113 114
115 1 8 2 114 115
116 1 8 2 115 116
117 1 8 2 116 117
118 1 8 2 117 118
119 1 8 2 118 119
120 1 8 2 119 120
121 1 8 2 120 121
122 1 8 2 121 122
123 1 8 2 122 123
124 1 8 2 123 124
125 1 8 2 124 125
126 1 8 2 125 126
127 1 8 2 126 6
128 1 10 2 7 8
129 1 11 2 5 7
130 1 12 2 6 8
1 15 1 1 1
2 15 2 1 2
3 15 15 1 7
4 15 16 1 8
$EndElements
#coding-Utf-8-*-
from gmshpy import *
#from dG3Dpy import *
from dG3DpyDebug import *
BulkLawNum1 = 1
ClLawNum1 = 2
lawnum2 = 3
lawnum3 = 4
# material law-----------------------------------------------
rho = 1400.0
E = 231.2e9
cl = 2.33e-5 * 2.33e-5
nu = 0.2
Cl_Law1 = IsotropicCLengthLaw(ClLawNum1,cl)
BulkLaw1 = NeoHookeanElasticPotential(BulkLawNum1, E, nu)
law1 = dG3DHyperelasticMaterialLaw(lawnum2,rho,BulkLaw1)
# BulkLaw1.setNonLocalLimitingMethod(0)
law1.setUseBarF(False)
mylaw1 = GenericCrackPhaseFieldDG3DMaterialLaw(lawnum3, rho, Cl_Law1, 90.0)
mylaw1.setMechanicalMaterialLaw(law1)
# domain
myfield1 = 120
fulldg = False
beta1 = 30.
eqRatio = 1.e8
space1 = 0 # function space (Lagrange=0)
mydom = dG3DDomain(10,myfield1,space1,lawnum3,fulldg,2,1)
# mydom.setNonLocalEqRatio(eqRatio)
mydom.stabilityParameters(30.)
PlaneStress = False #True
mydom.setPlaneStressState(PlaneStress)
#mydom.matrixByPerturbation(1,1,1,1e-8)
mydom.usePrimaryShapeFunction(3)
# solver
sol = 2 # Gmm=0 (default) Taucs=1 PETsc=2
soltype = 1 # StaticLinear=0 (default) StaticNonLinear=1
nstepImpl = 50 # number of step (used only if soltype=1)
nstepArch= 1 # Number of step between 2 archiving (used only if soltype=1)
ftime =5.e-2 # Final time (used only if soltype=1)
tol=1.0e-6
absTol=1.e-9
mysolver = nonLinearMechSolver(10)
geofile = 'Sample.geo'
meshfile = 'Sample.msh'
mysolver.createModel(geofile,meshfile,2,1)
mysolver.addMaterialLaw(mylaw1)
mysolver.addDomain(mydom)
mysolver.Scheme(soltype)
mysolver.Solver(sol)
mysolver.snlData(nstepImpl,ftime,tol,absTol)
mysolver.snlManageTimeStep(50, 3, 2, 10)
mysolver.stepBetweenArchiving(nstepArch)
# BC
mysolver.displacementBC("Edge",100,0,0.)
# mysolver.displacementBC("Face",1041,4,0.)
# mysolver.displacementBC("Face",1041,3,0.)
mysolver.displacementBC("Face",120,2,0.)
mysolver.displacementBC("Node",1011,1,0.)
mysolver.displacementBC("Node",1015,1,0.)
mysolver.displacementBC("Edge",110,1,0.)
mysolver.sameDisplacementBC("Edge",111,1016,1)
# method=0 # method 0
# mysolver.pathFollowing(True,method)
# if method==0:
# mysolver.setPathFollowingIncrementAdaptation(True,4,0.3)
# mysolver.setPathFollowingControlType(0)
# mysolver.setPathFollowingCorrectionMethod(0)
# mysolver.setPathFollowingArcLengthStep(1e-1)
# mysolver.setBoundsOfPathFollowingArcLengthSteps(0,0.1);
# elif method==1:
# mysolver.setPathFollowingIncrementAdaptationLocal(False,5,7)
# mysolver.setPathFollowingLocalSteps(10.,1e-3)
# mysolver.setPathFollowingLocalIncrementType(1);
# mysolver.setPathFollowingSwitchCriterion(0.3)
# mysolver.setBoundsOfPathFollowingLocalSteps(1e-13,4.)
# elif method==2:
# mysolver.setPathFollowingIncrementAdaptation(True,4,0.3)
# mysolver.setPathFollowingControlType(0)
# mysolver.setPathFollowingCorrectionMethod(0)
# mysolver.setPathFollowingArcLengthStep(1e-1)
# mysolver.setBoundsOfPathFollowingArcLengthSteps(0,2.);
# mysolver.sameDisplacementBC("Edge",101,1015,0)
# mysolver.forceBC("Node",1015,0,1.e4)
mysolver.displacementBC("Node",1015,0,0.5e-3)
# stopCri = EndSchemeMonitoringWithZeroForceOnPhysical(0.1,"Edge",101,0)
# mysolver.endSchemeMonitoring(stopCri)
# needed resultes
# mysolver.internalPointBuildView("svm",IPField.SVM, 1, 1)
mysolver.internalPointBuildView("sig_xx",IPField.SIG_XX, 1, 1)
mysolver.internalPointBuildView("sig_yy",IPField.SIG_YY, 1, 1)
#mysolver.internalPointBuildView("sig_zz",IPField.SIG_ZZ, 1, 1)
mysolver.internalPointBuildView("sig_xy",IPField.SIG_XY, 1, 1)
#mysolver.internalPointBuildView("sig_yz",IPField.SIG_YZ, 1, 1)
#mysolver.internalPointBuildView("sig_xz",IPField.SIG_XZ, 1, 1)
# mysolver.internalPointBuildView("epsilon_xx",IPField.STRAIN_XX, 1, 1)
# mysolver.internalPointBuildView("epsilon_yy",IPField.STRAIN_YY, 1, 1)
#mysolver.internalPointBuildView("epsilon_zz",IPField.STRAIN_ZZ, 1, 1)
# mysolver.internalPointBuildView("epsilon_xy",IPField.STRAIN_XY, 1, 1)
#mysolver.internalPointBuildView("epsilon_yz",IPField.STRAIN_YZ, 1, 1)
#mysolver.internalPointBuildView("epsilon_xz",IPField.STRAIN_XZ, 1, 1)
# mysolver.internalPointBuildView("epl",IPField.PLASTICSTRAIN, 1, 1)
mysolver.internalPointBuildView("damage",IPField.DAMAGE, 1, 1)
# mysolver.internalPointBuildView("FiberDamage",IPField.INC_DAMAGE, 1, 1)
# mysolver.archivingIPOnPhysicalGroup("Face",myfield1, IPField.DAMAGE,IPField.MIN_VALUE, 1);
# mysolver.archivingIPOnPhysicalGroup("Face",myfield1, IPField.DAMAGE,IPField.MAX_VALUE, 1);
# mysolver.archivingIPOnPhysicalGroup("Face",myfield1, IPField.DAMAGE,IPField.MEAN_VALUE, 1);
# mysolver.archivingIPOnPhysicalGroup("Face",myfield1, IPField.INC_DAMAGE,IPField.MIN_VALUE, 1);
# mysolver.archivingIPOnPhysicalGroup("Face",myfield1, IPField.INC_DAMAGE,IPField.MAX_VALUE, 1);
# mysolver.archivingIPOnPhysicalGroup("Face",myfield1, IPField.INC_DAMAGE,IPField.MEAN_VALUE, 1);
mysolver.archivingForceOnPhysicalGroup("Edge", 101, 0, 1)
mysolver.archivingForceOnPhysicalGroup("Edge", 100, 0, 1)
mysolver.archivingNodeDisplacement(1015,0,1)
mysolver.solve()
# check = TestCheck()
# check.equal(3.645579e-05,mysolver.getArchivedNodalValue(1015,0,mysolver.displacement),1.e-6)
// Plane strain specimen tests with notch
// Geometry parameters
//====================================================
// Units
mm=1.0e-3;
L = 1.*mm;
L_notch = L/2.;
tol = 1.*mm*mm;
LoadingMode = 0;
// Mesh size
//====================================================
sl1 = L/2.;
sl2 = L/100;
sl3 = L/10;
// Points
//====================================================
Point(1)= {0.,0.,0.,sl1};
Point(3)={L,L,0.,sl3};
Point(4)={0.,L,0.,sl1};
Point(5)={0.,L/2.+tol,0.,sl1};
Point(6)={L/2.,L/2.,0.,sl2};
Point(7)={0.,L/2.-tol,0.,sl1};
If (LoadingMode == 0)
Point(21)= {L*0.6,0.,0.,sl1};
Point(2)={L,0.,0.,sl3};
Point(22)= {L,L*0.43,0.,sl2};
Point(23)= {L,L*0.57,0.,sl2};
Else
Point(21)= {L*0.8,0.,0.,sl2};
Point(2)={L,0.,0.,sl2};
Point(22)= {L,L*0.6,0.,sl2*2.};
Point(23)= {L,L*0.8,0.,sl1};
EndIf
//
//
//
//// Lines
////====================================================
Line(1)={1,21};
Line(21) = {21,2};
Line(2)={2,22};
Line(22)={22,23};
Line(23)={23,3};
Line(3)={3,4};
Line(4)={4,5};
Line(5)={5,6};
Line(6)={6,7};
Line(7)={7,1};
// Surfaces
//====================================================
Line Loop(10) = {1,21,2,22,23,3,4,5,6,7};
Plane Surface(11) = {-10};
// Mesh size
//====================================================
If (LoadingMode == 0)
Point(30) = {L*0.5,L*0.43,0.,sl2};
Point(31) = {L*0.5,L*0.57,0.,sl2};
Point(32) = {L*0.6,L*0.43,0.,sl2};
Point(33) = {L*0.6,L*0.57,0.,sl2};
Point(34) = {L*0.7,L*0.43,0.,sl2};
Point(35) = {L*0.7,L*0.57,0.,sl2};
Point(36) = {L*0.81,L*0.43,0.,sl2};
Point(37) = {L*0.79,L*0.57,0.,sl2};
Point(38) = {L*0.89,L*0.43,0.,sl2};
Point(39) = {L*0.91,L*0.57,0.,sl2};
//Point(11) = {L*0.25,L*0.75,0.,sl2};
Point{30,31,32,33,34,35,36,37,38,39} In Surface{11};
Else
Point(29) = {L*0.45,L*0.48,0.,sl2};
Point(30) = {L*0.45,L*0.4,0.,sl2};
Point(31) = {L*0.5,L*0.3,0.,sl2};
Point(32) = {L*0.6,L*0.2,0.,sl2};
Point(33) = {L*0.7,L*0.1,0.,sl2};
Point(34) = {L*0.6,L*0.50,0.,sl2};
Point(35) = {L*0.7,L*0.35,0.,sl2};
Point(36) = {L*0.8,L*0.25,0.,sl2};
Point(37) = {L*0.9,L*0.15,0.,sl2};
Point(38) = {L*0.5,L*0.52,0.,sl2};
Point(39) = {L*0.9,L*0.55,0.,sl2*2.};
//Point(39) = {L*0.85,L*0.5,0.,sl2};
Point{29,30,31,32,33,34,35,36,37,38,39} In Surface{11};
EndIf
Extrude {0.0 , 0.0 , 0.25e-3} {Surface{11}; Layers{1}; Recombine;}
// Physical
//====================================================
Physical Volume(1000) = {1};
Physical Surface(100) = {11};
Physical Surface(200) = {54};
Physical Surface(300) = {58,62,66};
Physical Surface(400) = {70,74};
Physical Surface(500) = {75};
Physical Surface(600) = {38, 50};
Physical Line(101) = {1,21};
Physical Line(102) = {2,22};
Physical Line(103) = {3};
Physical Line(104) = {4,7};
Physical Point(111) = {3};
Physical Point(112) = {2};
This diff is collapsed.
// Plane strain specimen tests with notch
SetFactory('OpenCASCADE');
Mesh.SecondOrderLinear = 1;
// Geometry parameters
//====================================================
// Units
mm=1.0e-3;
L = 5.*mm;
L_notch = 2.*mm;
H = 10.*mm;
W = 2.*mm;
tol = 5.*mm*mm;
LoadingMode = 0;
// Mesh size
//====================================================
sl1 = H/2.;
sl2 = H/100;
sl3 = H/10;
// Points
//====================================================
Point(1)= {0.,0.,0.,sl1};
Point(3)={L,H,0.,sl3};
Point(4)={0.,H,0.,sl1};
Point(5)={0.,H/2.+tol,0.,sl1};
Point(6)={L_notch,H/2.,0.,sl2};
Point(7)={0.,H/2.-tol,0.,sl1};
Point(8)= {0.,0.,W,sl1};
Point(10)={L,H,W,sl3};
Point(11)={0.,H,W,sl1};
Point(12)={0.,H/2.+tol,W,sl1};
Point(13)={L_notch,H/2.,W,sl2};
Point(14)={0.,H/2.-tol,W,sl1};
If (LoadingMode == 0)
Point(21)= {L*0.6,0.,0.,sl1};
Point(2)={L,0.,0.,sl3};
Point(22)= {L,H*0.43,0.,sl2};
Point(23)= {L,H*0.57,0.,sl2};
Point(31)= {L*0.6,0.,W,sl1};
Point(9)={L,0.,W,sl3};
Point(32)= {L,H*0.43,W,sl2};
Point(33)= {L,H*0.57,W,sl2};
Else
Point(21)= {L*0.8,0.,0.,sl2};
Point(2)={L,0.,0.,sl2};
Point(22)= {L,L*0.6,0.,sl2*2.};
Point(23)= {L,L*0.8,0.,sl1};
EndIf
// Lines
//====================================================
Line(1)={1,21};
Line(21) = {21,2};
Line(2)={2,22};
Line(22)={22,23};
Line(23)={23,3};
Line(3)={3,4};
Line(4)={4,5};
Line(5)={5,6};
Line(6)={6,7};
Line(7)={7,1};
Line(8)={8,31};
Line(31) = {31,9};
Line(9)={9,32};
Line(32)={32,33};
Line(33)={33,10};
Line(10)={10,11};
Line(11)={11,12};
Line(12)={12,13};
Line(13)={13,14};
Line(14)={14,8};
Line(41)={1,8};
Line(42)={21,31};
Line(43)={2,9};
Line(44)={22,32};
Line(45)={23,33};
Line(46)={3,10};
Line(47)={4,11};
Line(48)={5,12};
Line(49)={6,13};
Line(50)={7,14};
//// Surfaces
////====================================================
//
Line Loop(100) = {1,21,2,22,23,3,4,5,6,7};
Plane Surface(11) = {-100};
Line Loop(200) = {8,31,9,32,33,10,11,12,13,14};
Plane Surface(22) = {-200};
Line Loop(300) = {41,8,-42,-1};
Plane Surface(33) = {-300};
Line Loop(400) = {42,31,-43,-21};
Plane Surface(44) = {-400};
Line Loop(500) = {43,9,-44,-2};
Plane Surface(55) = {-500};
Line Loop(600) = {44,32,-45,-22};
Plane Surface(66) = {-600};
Line Loop(700) = {45,33,-46,-23};
Plane Surface(77) = {-700};
Line Loop(800) = {47,-10,-46,3};
Plane Surface(88) = {-800};
Line Loop(900) = {48,-11,-47,4};
Plane Surface(99) = {-900};
Line Loop(910) = {48,12,-49,-5};
Plane Surface(100) = {-910};
Line Loop(920) = {50,-13,-49,6};
Plane Surface(101) = {-920};
Coherence;
Line Loop(930) = {41,-14,-50,7};
Plane Surface(105) = {-930};
Coherence;
Surface Loop(500) = {33, 44, 55, 66, 77, 88, 99, 100, 101, 105, 11, 22};
Volume(20) = {500};
Coherence;
//
//
//// Mesh size
////====================================================
//
//If (LoadingMode == 0)
// Point(30) = {L*0.5,L*0.43,0.,sl2};
// Point(31) = {L*0.5,L*0.57,0.,sl2};
// Point(32) = {L*0.6,L*0.43,0.,sl2};
// Point(33) = {L*0.6,L*0.57,0.,sl2};
// Point(34) = {L*0.7,L*0.43,0.,sl2};
// Point(35) = {L*0.7,L*0.57,0.,sl2};
// Point(36) = {L*0.81,L*0.43,0.,sl2};
// Point(37) = {L*0.79,L*0.57,0.,sl2};
// Point(38) = {L*0.89,L*0.43,0.,sl2};
// Point(39) = {L*0.91,L*0.57,0.,sl2};
// //Point(11) = {L*0.25,L*0.75,0.,sl2};
// Point{30,31,32,33,34,35,36,37,38,39} In Surface{11};
//
//
//Else
// Point(29) = {L*0.45,L*0.48,0.,sl2};
//
// Point(30) = {L*0.45,L*0.4,0.,sl2};
// Point(31) = {L*0.5,L*0.3,0.,sl2};
// Point(32) = {L*0.6,L*0.2,0.,sl2};
// Point(33) = {L*0.7,L*0.1,0.,sl2};
//
// Point(34) = {L*0.6,L*0.50,0.,sl2};
// Point(35) = {L*0.7,L*0.35,0.,sl2};
// Point(36) = {L*0.8,L*0.25,0.,sl2};
// Point(37) = {L*0.9,L*0.15,0.,sl2};
//
// Point(38) = {L*0.5,L*0.52,0.,sl2};
//
// Point(39) = {L*0.9,L*0.55,0.,sl2*2.};
// //Point(39) = {L*0.85,L*0.5,0.,sl2};
// Point{29,30,31,32,33,34,35,36,37,38,39} In Surface{11};
//EndIf
//
//
//Extrude {0.0 , 0.0 , 0.25e-3} {Surface{11}; Layers{1}; Recombine;}
//
//// Physical
////====================================================
//
Physical Volume(1000) = {20};
//
Physical Surface(100) = {11};
Physical Surface(200) = {88};
Physical Surface(300) = {55,66,77};
Physical Surface(400) = {33,44};
Physical Surface(500) = {22};
Physical Surface(600) = {99, 105};
//Physical Line(101) = {1,21};
//Physical Line(102) = {2,22};
//Physical Line(103) = {3};
//Physical Line(104) = {4,7};
//
//Physical Point(111) = {3};
//Physical Point(112) = {2};
This diff is collapsed.
#coding-Utf-8-*-
from gmshpy import*
from dG3DpyDebug import*
#from dG3Dpy import*
# Script for plane plate problem
# Options
# ===================================================================================
Flag_Fracture = 1 # Activate cracks and fracture laws ? (1=yes, 0=no)
Solving_mode = 1 # QStatic (=1) vs. Implicit (=2)
LoadingMode = 0
# Material law creation
# ===================================================================================
# numbering of material law (has to be unique for each law !)
BulkLawNum1 = 1
ClLawNum1 = 2
DamLawNum1 = 3
HardenLawNum1 = 4
InterfaceLawNum1 = 5
FracLawNum1 = 6
lawnum2 = 7
lawnum3 = 8
# material parameters - bulk law
rho = 2450.0 # kg/m3 : density
E = 210.0e9 # Pa : Young modulus
nu = 0.3 # - : poisson coefficient
cl = (0.015*1.0e-3) * (0.015*1.0e-3) # m2: non-local parameter (== NL-length^2)
# # material parameters (2) - cohesive law
# sigmac = 0.028*E #0.09*E # Pa critical effective stress value
# Dc = 0.9 # - critical damage value
# damageCheck = bool(0) # bool actived crack insertion by damage value
# h = 5.4*0.015*1.10*1e-3 # m band thickness
# mu = -10000. # - friction coefficient
# fmin = 1.00 # - spread min bond
# fmax =1.01 # - spread max bond
# betac = 0.87 # -
# Kp = 5.e12 # Pa/m2 penatly parameter (F = -K_p u_n^2)
# # ==O(E * beta_DG / mesh_size^2)
# GradFactor = 0.5 # - multiply grad of jump
# useFc = bool(1) # - activate normal comp. blocking of F
# delta0 = 1.0e-6
# material law creation
Cl_Law1 = IsotropicCLengthLaw(ClLawNum1,cl)
# Dam_Law1 = PowerDamageLaw(DamLawNum1, pi, pc, alpha, beta)
# #Dam_Law1 =ExponentialDamageLaw(DamLawNum1,1.e-4,0.1,500.)
# BulkLaw1 = NonLocalDamageIsotropicElasticityDG3DMaterialLaw(BulkLawNum1,rho,E,nu,Cl_Law1,Dam_Law1)
BulkLaw1 = biLogarithmicElasticPotential(BulkLawNum1, E, nu, 1)
# BulkLaw1 = NeoHookeanElasticPotential(BulkLawNum1, E, nu)
law1 = dG3DHyperelasticMaterialLaw(lawnum2,rho,BulkLaw1)
# BulkLaw1.setNonLocalLimitingMethod(0)
law1.setUseBarF(False)
# InterfaceLaw1 = CohesiveBand3DLaw(InterfaceLawNum1,sigmac,Dc,h,mu,fmin,fmax, betac,Kp,damageCheck,GradFactor,useFc,delta0)
#InterfaceLaw1 = NonLocalDamageLinearCohesive3DLaw(InterfaceLawNum1,21.4*1.0e3,sigmac,Dc,betac,fmin, fmax,1e7,damageCheck);
pi = 0.01
b = 0.
K = 20.
# FracLaw1 = FractureByCohesive3DLaw(FracLawNum1,BulkLawNum1,InterfaceLawNum1)
damlaw = ExponentialDamageLaw(1, pi, b, K)
mylaw1 = NonLocalDamageHyperelasticDG3DMaterialLaw(lawnum3, rho, BulkLaw1, Cl_Law1, damlaw)
# mylaw1 = GenericCrackPhaseFieldDG3DMaterialLaw(lawnum3, rho, Cl_Law1, 270.0, 1.0)
# mylaw1.setMechanicalMaterialLaw(law1)
# Solver parameters
# ===================================================================================
sol = 2 # Library for solving: Gmm=0 (default) Taucs=1 PETsc=2
if Solving_mode == 1:
soltype = 1 # StaticLinear=0 (default, StaticNonLinear=1, Explicit=2,
# Multi=3, Implicit=4, Eigen=5)
nstep = 100 # Number of step
ftime =1.0 # Final time
tol=1.e-5 # Relative tolerance for NR scheme
tolAbs = 1.e-6 # Absolute tolerance for NR scheme
nstepArch= 1 # Number of step between 2 archiving
nstepArchEnergy = 1 # Number of step between 2 energy computation and archiving
nstepArchForce = 1 # Number of step between 2 force archiving
MaxIter = 200 # Maximum number of iterations
elif Solving_mode == 2:
# soltype = 3 # StaticLinear=0 (default, StaticNonLinear=1, Explicit=2,
# Multi=3, Implicit=4, Eigen=5
soltype = 4 #@Mohib changed to 4 from 2
nstep = 500 # Number of step between implicit iterations
nstepExpl = 2500 # Number of step between time step evaluation
# ftime =(1.0e-4)*0.1 # Final time
ftime =1.0e-2 # Final time
tol=1.e-4 # Relative tolerance for NR scheme
tolAbs = 1.e-6 # Absolute tolerance for NR scheme
nstepArch= 1 # Number of step between 2 archiving
nstepArchEnergy = 1 # Number of step between 2 energy computation and archiving
nstepArchForce = 1 # Number of step between 2 force archiving
MaxIter = 200 # Maximum number of iterations
StepIncrease = 2 # Number of successfull timestep before reincreasing it
StepReducFactor = 5.0 # Timestep reduction factor
NumberReduction = 4 # Maximum number of timespep reduction: max reduction = pow(StepReducFactor,this)
fullDg = False # O = CG, 1 = DG
dgnl = False # DG for non-local variables inside a domain (only if fullDg)
eqRatio = 1.0 #1.0e6 @Mohib changed from a large number # Ratio between "elastic" and non-local equations
space1 = 0 # Function space (Lagrange=0)
beta1 = 30.0 # Penality parameter for DG
# Domain creation
## ===================================================================================
# Example : field1 = dG3DDomain(numDomain, numPhysVolume, space, LawNum, fullDgFlag, dim=3, numNonlocalVar=1)
# with: - numDomain : domain tag
# - numPhysVolume : number of physical entities
# - space : function space (Lagrange=0)
# - LawNum : number of material law
# - fullDgFlag : O = CG, 1 = DG
# - nonLocalEqRatio : ratio between elastic and non-local system
# - dim : dimension of displacment
# - numNonlocalVar : numbver of non-local variable
numPhysVolume1 = 1000 # Number of a physical volume of the model in .geo
numDomain1 = 1001 # Number of the domain
# if Flag_Fracture == 0 :
# field1 = dG3DDomain(numDomain1,numPhysVolume1,space1,BulkLawNum1,fullDg,2,1)
# elif
# Flag_Fracture == 1 :
# field1 = dG3DDomain(numDomain1,numPhysVolume1,space1,FracLawNum1,fullDg,2,1)
field1 = dG3DDomain(numDomain1,numPhysVolume1,space1,lawnum3,fullDg,3,1)
field1.stabilityParameters(beta1) # Adding stability parameters (for DG)
field1.setNonLocalStabilityParameters(beta1,dgnl) # Adding stability parameters (for DG)
field1.setNonLocalEqRatio(eqRatio)
# field1.gaussIntegration(0,-1,-1)
# field1.setBulkDamageBlockedMethod(1)
# field1.forceCohesiveInsertionAllIPs(bool(1),0.7)
field1.usePrimaryShapeFunction(3)
# Solver creation
# ===================================================================================
mysolver = nonLinearMechSolver(numDomain1) # Solver associated with numSolver1
# Geometry and mesh loading
# ===================================================================================
meshfile="SEN3.msh" # name of mesh file
mysolver.loadModel(meshfile) # add mesh
mysolver.addDomain(field1) # add domain
# mysolver.addMaterialLaw(BulkLaw1) # add material law
# mysolver.addMaterialLaw(InterfaceLaw1)
# mysolver.addMaterialLaw(FracLaw1)
mysolver.addMaterialLaw(mylaw1)
# solver parametrisation
if Solving_mode == 1 :
mysolver.Scheme(soltype) # solver scheme
mysolver.Solver(sol) # solver choice
mysolver.snlData(nstep,ftime,tol,tolAbs) # solver parameters
mysolver.snlManageTimeStep(MaxIter,StepIncrease,StepReducFactor,NumberReduction) #timestep
mysolver.lineSearch(bool(0)) # lineSearch activation
# solver archiving
mysolver.stepBetweenArchiving(nstepArch) # archiving frequency for IPField
# mysolver.energyComputation(nstepArchEnergy) # archiving frequency for energy
elif Solving_mode == 2 :
mysolver.Scheme(soltype) # solver scheme
mysolver.Solver(sol) # solver choice
mysolver.snlData(nstep,ftime,tol,tolAbs) # solver parameters
mysolver.snlManageTimeStep(MaxIter,StepIncrease,StepReducFactor,NumberReduction) #timestep
mysolver.lineSearch(bool(0)) # LineSearch activation
# mysolver.explicitSpectralRadius(ftime,0.5,0.0) # Explicit integration parameters
mysolver.implicitSpectralRadius(0.05) # Implicit
# mysolver.explicitTimeStepEvaluation(nstepExpl)
# mysolver.addSystem(3,2) # adding system (numb. of var., solv. scheme)
# mysolver.addSystem(1,1)
# solver archiving
mysolver.stepBetweenArchiving(nstepArch) # archiving frequency
mysolver.energyComputation(nstepArchEnergy) # archiving frequency for energy
stiffModif = BFGSStiffnessModificationMonitoring(20, False) # stiffness will be re-computed after 20 iterations in each step
mysolver.stiffnessModification(stiffModif)
# Boundary conditions
# ===============================
# Example with constrained displacement : mysolver.displacementBC("type", numphys, comp, value)
# with: - type : Node / Edge / Face / Volume
# - numphys : number of physical entities
# - comp : number of variables on which it's applied (0 = u_x; 1 = u_y; 2 = u_z)
# - value : value at the end of the simulation
# NB: if explicit solving: displacment is replaced by a speed !
# Example with applied forces : mysolver.forceBC("type", numphys, comp, value)
if LoadingMode == 0:
#tot_disp = 0.01* 1.0*1.e-3
mysolver.displacementBC("Face",100,2,0.)
mysolver.displacementBC("Face",500,2,0.)
mysolver.displacementBC("Face",300,0,0.)
mysolver.displacementBC("Face",400,1,0.)
# face x = 0
# mysolver.displacementBC("Edge",101,1,0.) # face x = 0
# mysolver.displacementBC("Face",100,2,0.)
mysolver.forceBC("Face",200,1, 1.0e7)
mysolver.constraintBC("Face",200,1)
# #mysolver.displacementBC("Edge",103,1,tot_disp/ftime) # face x = L
# mysolver.displacementBC("Edge",103,0,0.) # face x = L
#mysolver.displacementBC("Edge",104,0,0.) # face y = 0
method = 2
mysolver.pathFollowing(True,method)
if method==0:
mysolver.setPathFollowingIncrementAdaptation(True,4,0.3)
mysolver.setPathFollowingControlType(0)
mysolver.setPathFollowingCorrectionMethod(0)
mysolver.setPathFollowingArcLengthStep(1e-1)
mysolver.setBoundsOfPathFollowingArcLengthSteps(0,0.1);
elif method==1:
mysolver.setPathFollowingIncrementAdaptationLocal(False,5,7)
mysolver.setPathFollowingLocalSteps(10.,1e-3)
mysolver.setPathFollowingLocalIncrementType(1);
mysolver.setPathFollowingSwitchCriterion(0.3)
mysolver.setBoundsOfPathFollowingLocalSteps(1e-13,4.)
elif method==2:
mysolver.setPathFollowingIncrementAdaptation(True,4,0.3)
mysolver.setPathFollowingControlType(0)
mysolver.setPathFollowingCorrectionMethod(0)
mysolver.setPathFollowingArcLengthStep(1e-1)
mysolver.setBoundsOfPathFollowingArcLengthSteps(0,2.);
# mysolver.pathFollowing(True,1)
# mysolver.setPathFollowingIncrementAdaptation(False,4)
# mysolver.setPathFollowingLocalSteps(9.0e-1,1.0e-3)
# mysolver.setPathFollowingLocalIncrementType(1);
elif LoadingMode == 1:
tot_disp = 0.02* 1.0*1.e-3
mysolver.displacementBC("Edge",101,0,0.) # face x = 0
mysolver.displacementBC("Edge",101,1,0.) # face x = 0
mysolver.displacementBC("Face",100,2,0.)
#mysolver.forceBC("Edge",103,0, 1.e8)
#mysolver.constraintBC("Edge",103,0)
mysolver.displacementBC("Edge",103,0,tot_disp/ftime) # face x = L
mysolver.displacementBC("Edge",103,1,0.) # face x = L
mysolver.displacementBC("Edge",102,1,0.)
mysolver.displacementBC("Edge",104,1,0.)
elif LoadingMode == 2:
tot_disp = 7.0e-3 * 1.e-3
mysolver.displacementBC("Face",100,2,0.) # face x = 0
mysolver.displacementBC("Face",500,2,0.) # face x = 0
mysolver.displacementBC("Face",300,0,0.) # face x = 0
mysolver.displacementBC("Face",400,1,0.) # face x = 0
mysolver.displacementBC("Face",200,1, tot_disp/ftime)
mysolver.displacementBC("Face",200,0, 0.)
# Variable storage
# ===============================
# Example : mysolver.internalPointBuildView("file_name", "comp", nbstep, ev)
# with: - file_name : the name of variable storage file
# - comp : saved values (IPField.SVM, IPField.SIG_XX, IPField::STRAIN_XZ,...
# see ipField.h for possibilities)
# - nbstep : storage frequency ?
# - ev : value type (default : mean)
#
# Example : mysolver.archivingForceOnPhysicalGroup("type", numphys, comp, nstep)
# with: - type : Node / Edge / Face / Volume
# - numphys : number of physical entities
# - comp : number of variables on which it's applied (0 = u_x; 1 = u_y; 2 = u_z)
# - nstep : step number between force archivingmysolver.internalPointBuildView("svm",IPField.SVM, 1, 1)
#
mysolver.internalPointBuildView("svm",IPField.SVM, 1, 1)
mysolver.internalPointBuildView("sig_xx",IPField.SIG_XX, 1, 1)
mysolver.internalPointBuildView("sig_yy",IPField.SIG_YY, 1, 1)
mysolver.internalPointBuildView("sig_zz",IPField.SIG_ZZ, 1, 1)
mysolver.internalPointBuildView("sig_xy",IPField.SIG_XY, 1, 1)
mysolver.internalPointBuildView("sig_yz",IPField.SIG_YZ, 1, 1)
mysolver.internalPointBuildView("sig_xz",IPField.SIG_XZ, 1, 1)
mysolver.internalPointBuildView("damage",IPField.DAMAGE,1,1)
# mysolver.internalPointBuildView("eff_eps",IPField.LOCAL_EQUIVALENT_STRAINS,1,1)
# mysolver.internalPointBuildView("NL_eps",IPField.NONLOCAL_EQUIVALENT_STRAINS,1,1)
# mysolver.internalPointBuildView("Blocked_element",IPField.DAMAGE_IS_BLOCKED,1,1)
mysolver.archivingForceOnPhysicalGroup("Edge", 103, 0, nstepArchForce)
mysolver.archivingForceOnPhysicalGroup("Edge", 103, 1, nstepArchForce)
mysolver.archivingForceOnPhysicalGroup("Edge", 101, 1, nstepArchForce)
mysolver.archivingForceOnPhysicalGroup("Edge", 101, 0, nstepArchForce)
mysolver.archivingNodeDisplacement(111,0, nstepArchForce)
mysolver.archivingNodeDisplacement(111,1, nstepArchForce)
# Solving
# ===========
mysolver.solve()
//Fiber reinforced composite cross-section
Geometry.AutoCoherence= 1.0e-9;
//defination of unit
unit = 1.0e-3;
//characteristic mesh size
h_x=5*unit;
h_y=5*unit;//0.1
h_z = 0.02*unit;
L = 0.7*unit;
nx= 1;
ny= 1;
cl=5.0*unit;
//point
Point(1) = {0.0 , 0.0 , 0.0, cl};
Point(2) = {0.0 , ny*h_y , 0.0, cl};
Line(1) = {1,2};
Extrude {nx*h_x , 0.0 , 0.0} {Line{1}; Layers{nx}; Recombine;}
//Extrude {L*2 , 0.0 , 0.0} {Line{2};Layers{60}; Recombine;}
//Extrude {nx*h_x , 0.0 , 0.0} {Line{6}; Layers{nx}; Recombine;}
//
//
Physical Line(101) = {3};
Physical Line(102) = {2};
Physical Line(103) = {4};
Physical Line(104) = {1};
Physical Surface(120) = {5};
//
Physical Point(201) = {1};
Physical Point(202) = {3};
Physical Point(203) = {4};
Physical Point(204) = {2};
//Physical Point(1016) = {16};
Coherence;
$MeshFormat
4.1 0 8
$EndMeshFormat
$Entities
4 4 1 0
1 0 0 0 1 201
2 0 0.005 0 1 204
3 0.005 0 0 1 202
4 0.005 0.005 0 1 203
1 0 0 0 0 0.005 0 1 104 2 1 -2
2 0.005 0 0 0.005 0.005 0 1 102 2 3 -4
3 0 0 0 0.005 0 0 1 101 2 1 -3
4 0 0.005 0 0.005 0.005 0 1 103 2 2 -4
5 0 0 0 0.005 0.005 0 1 120 4 1 4 -2 -3
$EndEntities
$Nodes
9 9 1 9
0 1 0 1
1
0 0 0
0 2 0 1
2
0 0.005 0
0 3 0 1
3
0.005 0 0
0 4 0 1
4
0.005 0.005 0
1 1 0 1
5
0 0.002499999999994851 0
1 2 0 1
6
0.005 0.002499999999994851 0
1 3 0 1
7
0.002499999999994851 0 0
1 4 0 1
8
0.002499999999994851 0.005 0
2 5 0 1
9
0.002499999999994851 0.002499999999994851 0
$EndNodes
$Elements
9 9 1 9
0 1 15 1
1 1
0 2 15 1
2 2
0 3 15 1
3 3
0 4 15 1
4 4
1 1 8 1
5 1 2 5
1 2 8 1
6 3 4 6
1 3 8 1
7 1 3 7
1 4 8 1
8 2 4 8
2 5 10 1
9 1 2 4 3 5 8 6 7 9
$EndElements
#coding-Utf-8-*-
from gmshpy import *
from dG3DpyDebug import*
#script to launch beam problem with a python script
# material law
# lawnum = 1 # unique number of law
bulk_law_num = 1
mech_law_num = 2
cl_law_num = 3
pf_law_num = 4
rho = 7850
# young = 28.9e9
E = 28.9e9
nu = 0.3
cl_ = 5.0e-3 * 5.0e-3 * 1.1 * 1.1
# geometry
meshfile="Sample.msh" # name of mesh file
# solver
sol = 2 # Gmm=0 (default) Taucs=1 PETsc=2
soltype = 1 # StaticLinear=0 (default) StaticNonLinear=1
nstep = 100 # number of step (used only if soltype=1)
ftime =1. # Final time (used only if soltype=1)
tol=1.e-6 # relative tolerance for NR scheme (used only if soltype=1)
nstepArch=1 # Number of step between 2 archiving (used only if soltype=1)
fullDg = 0 #O = CG, 1 = DG
space1 = 0 # function space (Lagrange=0)
beta1 = 100
# compute solution and BC (given directly to the solver
#elasticLaw = smallStrainLinearElasticPotential(1,young,nu)
bulk_law = NeoHookeanElasticPotential(bulk_law_num, E, nu)
#mech_law = dG3DHyperelasticMaterialLaw(mech_law_num, rho, bulk_law)
#mech_law.setUseBarF(False)
# elasticLaw = biLogarithmicElasticPotential(1,young,nu,1)
cl = IsotropicCLengthLaw(cl_law_num, cl_)
pi = 0.01
b = 0.
K = 20.
damlaw = ExponentialDamageLaw(1, pi, b, K)
pf_law = NonLocalDamageHyperelasticDG3DMaterialLaw(pf_law_num, rho, bulk_law, cl, damlaw)
# cohesive law
# creation of ElasticField
# nfield = 11 # number of the field (physical number of surface)
nfield = 120 # number of the field (physical number of surface)
myfield1 = dG3DDomain(1000,nfield,space1,pf_law_num,fullDg,2,1)
myfield1.stabilityParameters(beta1)
myfield1.setNonLocalStabilityParameters(beta1,True)
myfield1.setNonLocalEqRatio(1.e6)
myfield1.averageStrainBased(False)
PlaneStress = False #True
myfield1.setPlaneStressState(PlaneStress)
# creation of Solver
mysolver = nonLinearMechSolver(1000)
mysolver.loadModel(meshfile)
mysolver.addDomain(myfield1)
mysolver.addMaterialLaw(pf_law)
mysolver.Scheme(soltype)
mysolver.Solver(sol)
mysolver.snlData(nstep,ftime,tol)
mysolver.stepBetweenArchiving(nstepArch)
# BC
# mysolver.displacementBC("Volume",11,2,0.)
mysolver.displacementBC("Face",120,2,0.)
# mysolver.displacementBC("Face",12,1,0.)
mysolver.displacementBC("Edge",101,1,0.)
# mysolver.displacementBC("Face",14,0,0.)
mysolver.displacementBC("Edge",104,0,0.)
# mysolver.displacementBC("Face",13,1,1e-3)
mysolver.displacementBC("Edge",103,1,1e-3)
mysolver.internalPointBuildView("svm",IPField.SVM, 1, 1)
mysolver.internalPointBuildView("sig_xx",IPField.SIG_XX, 1, 1)
mysolver.internalPointBuildView("sig_yy",IPField.SIG_YY, 1, 1)
mysolver.internalPointBuildView("sig_zz",IPField.SIG_ZZ, 1, 1)
mysolver.internalPointBuildView("sig_xy",IPField.SIG_XY, 1, 1)
mysolver.internalPointBuildView("sig_yz",IPField.SIG_YZ, 1, 1)
mysolver.internalPointBuildView("sig_xz",IPField.SIG_XZ, 1, 1)
mysolver.internalPointBuildView("Defo Energy",IPField.DEFO_ENERGY, 1, 1)
mysolver.internalPointBuildView("Damage",IPField.DAMAGE, 1, 1)
# mysolver.archivingForceOnPhysicalGroup("Face", 12, 1)
mysolver.archivingForceOnPhysicalGroup("Edge", 101, 1, 1)
# mysolver.archivingNodeDisplacement(19,1,1)
mysolver.archivingNodeDisplacement(204,1,1)
mysolver.solve()
# check = TestCheck()
# check.equal(-1.294255e+02,mysolver.getArchivedForceOnPhysicalGroup("Face", 12, 1),1.e-4)
// Plane strain specimen tests with notch
// Geometry parameters
//====================================================
// Units
mm=1.0e-3;
L = 1.*mm;
L_notch = L/2.;
tol = 1.*mm*mm;
LoadingMode = 0;
// Mesh size
//====================================================
sl1 = L/2.;
sl2 = L/100;
sl3 = L/10;
// Points
//====================================================
Point(1)= {0.,0.,0.,sl1};
Point(3)={L,L,0.,sl3};
Point(4)={0.,L,0.,sl1};
Point(5)={0.,L/2.+tol,0.,sl1};
Point(6)={L/2.,L/2.,0.,sl2};
Point(7)={0.,L/2.-tol,0.,sl1};
If (LoadingMode == 0)
Point(21)= {L*0.6,0.,0.,sl1};
Point(2)={L,0.,0.,sl3};
Point(22)= {L,L*0.43,0.,sl2};
Point(23)= {L,L*0.57,0.,sl2};
Else
Point(21)= {L*0.8,0.,0.,sl2};
Point(2)={L,0.,0.,sl2};
Point(22)= {L,L*0.6,0.,sl2*2.};
Point(23)= {L,L*0.8,0.,sl1};
EndIf
// Lines
//====================================================
Line(1)={1,21};
Line(21) = {21,2};
Line(2)={2,22};
Line(22)={22,23};
Line(23)={23,3};
Line(3)={3,4};
Line(4)={4,5};
Line(5)={5,6};
Line(6)={6,7};
Line(7)={7,1};
// Surfaces
//====================================================
Line Loop(10) = {1,21,2,22,23,3,4,5,6,7};
Plane Surface(11) = {-10};
// Mesh size
//====================================================
If (LoadingMode == 0)
Point(30) = {L*0.5,L*0.43,0.,sl2};
Point(31) = {L*0.5,L*0.57,0.,sl2};
Point(32) = {L*0.6,L*0.43,0.,sl2};
Point(33) = {L*0.6,L*0.57,0.,sl2};
Point(34) = {L*0.7,L*0.43,0.,sl2};
Point(35) = {L*0.7,L*0.57,0.,sl2};
Point(36) = {L*0.81,L*0.43,0.,sl2};
Point(37) = {L*0.79,L*0.57,0.,sl2};
Point(38) = {L*0.89,L*0.43,0.,sl2};
Point(39) = {L*0.91,L*0.57,0.,sl2};
//Point(11) = {L*0.25,L*0.75,0.,sl2};
Point{30,31,32,33,34,35,36,37,38,39} In Surface{11};
Else
Point(29) = {L*0.45,L*0.48,0.,sl2};
Point(30) = {L*0.45,L*0.4,0.,sl2};
Point(31) = {L*0.5,L*0.3,0.,sl2};
Point(32) = {L*0.6,L*0.2,0.,sl2};
Point(33) = {L*0.7,L*0.1,0.,sl2};
Point(34) = {L*0.6,L*0.50,0.,sl2};
Point(35) = {L*0.7,L*0.35,0.,sl2};
Point(36) = {L*0.8,L*0.25,0.,sl2};
Point(37) = {L*0.9,L*0.15,0.,sl2};
Point(38) = {L*0.5,L*0.52,0.,sl2};
Point(39) = {L*0.9,L*0.55,0.,sl2*2.};
//Point(39) = {L*0.85,L*0.5,0.,sl2};
Point{29,30,31,32,33,34,35,36,37,38,39} In Surface{11};
EndIf
// Physical
//====================================================
Physical Surface(100) = {11};
Physical Line(101) = {1,21};
Physical Line(102) = {2,22};
Physical Line(103) = {3};
Physical Line(104) = {4,7};
Physical Point(111) = {3};
Physical Point(112) = {2};
This diff is collapsed.
#coding-Utf-8-*-
from gmshpy import*
from dG3DpyDebug import*
#from dG3Dpy import*
# Script for plane plate problem
# Options
# ===================================================================================
Flag_Fracture = 1 # Activate cracks and fracture laws ? (1=yes, 0=no)
Solving_mode = 1 # QStatic (=1) vs. Implicit (=2)
LoadingMode = 0
# Material law creation
# ===================================================================================
# numbering of material law (has to be unique for each law !)
BulkLawNum1 = 1
ClLawNum1 = 2
DamLawNum1 = 3
HardenLawNum1 = 4
InterfaceLawNum1 = 5
FracLawNum1 = 6
lawnum2 = 7
lawnum3 = 8
# material parameters - bulk law
rho = 2450.0 # kg/m3 : density
E = 210.0e9 # Pa : Young modulus
nu = 0.3 # - : poisson coefficient
cl = (0.015*1.0e-3)*(0.015*1.0e-3)*1.10*1.10 # m2: non-local parameter (== NL-length^2)
# # material parameters (2) - cohesive law
# sigmac = 0.028*E #0.09*E # Pa critical effective stress value
# Dc = 0.9 # - critical damage value
# damageCheck = bool(0) # bool actived crack insertion by damage value
# h = 5.4*0.015*1.10*1e-3 # m band thickness
# mu = -10000. # - friction coefficient
# fmin = 1.00 # - spread min bond
# fmax =1.01 # - spread max bond
# betac = 0.87 # -
# Kp = 5.e12 # Pa/m2 penatly parameter (F = -K_p u_n^2)
# # ==O(E * beta_DG / mesh_size^2)
# GradFactor = 0.5 # - multiply grad of jump
# useFc = bool(1) # - activate normal comp. blocking of F
# delta0 = 1.0e-6
# material law creation
Cl_Law1 = IsotropicCLengthLaw(ClLawNum1,cl)
# Dam_Law1 = PowerDamageLaw(DamLawNum1, pi, pc, alpha, beta)
# #Dam_Law1 =ExponentialDamageLaw(DamLawNum1,1.e-4,0.1,500.)
# BulkLaw1 = NonLocalDamageIsotropicElasticityDG3DMaterialLaw(BulkLawNum1,rho,E,nu,Cl_Law1,Dam_Law1)
BulkLaw1 = NeoHookeanElasticPotential(BulkLawNum1, E, nu)
law1 = dG3DHyperelasticMaterialLaw(lawnum2,rho,BulkLaw1)
# BulkLaw1.setNonLocalLimitingMethod(0)
law1.setUseBarF(False)
# InterfaceLaw1 = CohesiveBand3DLaw(InterfaceLawNum1,sigmac,Dc,h,mu,fmin,fmax, betac,Kp,damageCheck,GradFactor,useFc,delta0)
#InterfaceLaw1 = NonLocalDamageLinearCohesive3DLaw(InterfaceLawNum1,21.4*1.0e3,sigmac,Dc,betac,fmin, fmax,1e7,damageCheck);
# FracLaw1 = FractureByCohesive3DLaw(FracLawNum1,BulkLawNum1,InterfaceLawNum1)
mylaw1 = GenericCrackPhaseFieldDG3DMaterialLaw(lawnum3, rho, Cl_Law1, 2700.0, 1.0)
mylaw1.setMechanicalMaterialLaw(law1)
# Solver parameters
# ===================================================================================
sol = 2 # Library for solving: Gmm=0 (default) Taucs=1 PETsc=2
if Solving_mode == 1:
soltype = 1 # StaticLinear=0 (default, StaticNonLinear=1, Explicit=2,
# Multi=3, Implicit=4, Eigen=5)
nstep = 100 # Number of step
ftime =1.0 # Final time
tol=1.e-5 # Relative tolerance for NR scheme
tolAbs = 1.e-6 # Absolute tolerance for NR scheme
nstepArch= 1 # Number of step between 2 archiving
nstepArchEnergy = 1 # Number of step between 2 energy computation and archiving
nstepArchForce = 1 # Number of step between 2 force archiving
MaxIter = 200 # Maximum number of iterations
elif Solving_mode == 2:
# soltype = 3 # StaticLinear=0 (default, StaticNonLinear=1, Explicit=2,
# Multi=3, Implicit=4, Eigen=5
soltype = 4 #@Mohib changed to 4 from 2
nstep = 500 # Number of step between implicit iterations
nstepExpl = 2500 # Number of step between time step evaluation
# ftime =(1.0e-4)*0.1 # Final time
ftime =1.0e-2 # Final time
tol=1.e-4 # Relative tolerance for NR scheme
tolAbs = 1.e-6 # Absolute tolerance for NR scheme
nstepArch= 1 # Number of step between 2 archiving
nstepArchEnergy = 1 # Number of step between 2 energy computation and archiving
nstepArchForce = 1 # Number of step between 2 force archiving
MaxIter = 200 # Maximum number of iterations
StepIncrease = 2 # Number of successfull timestep before reincreasing it
StepReducFactor = 5.0 # Timestep reduction factor
NumberReduction = 4 # Maximum number of timespep reduction: max reduction = pow(StepReducFactor,this)
fullDg = False # O = CG, 1 = DG
dgnl = False # DG for non-local variables inside a domain (only if fullDg)
eqRatio = 1.0 #1.0e6 @Mohib changed from a large number # Ratio between "elastic" and non-local equations
space1 = 0 # Function space (Lagrange=0)
beta1 = 30.0 # Penality parameter for DG
# Domain creation
## ===================================================================================
# Example : field1 = dG3DDomain(numDomain, numPhysVolume, space, LawNum, fullDgFlag, dim=3, numNonlocalVar=1)
# with: - numDomain : domain tag
# - numPhysVolume : number of physical entities
# - space : function space (Lagrange=0)
# - LawNum : number of material law
# - fullDgFlag : O = CG, 1 = DG
# - nonLocalEqRatio : ratio between elastic and non-local system
# - dim : dimension of displacment
# - numNonlocalVar : numbver of non-local variable
numPhysVolume1 = 100 # Number of a physical volume of the model in .geo
numDomain1 = 1001 # Number of the domain
# if Flag_Fracture == 0 :
# field1 = dG3DDomain(numDomain1,numPhysVolume1,space1,BulkLawNum1,fullDg,2,1)
# elif
# Flag_Fracture == 1 :
# field1 = dG3DDomain(numDomain1,numPhysVolume1,space1,FracLawNum1,fullDg,2,1)
field1 = dG3DDomain(numDomain1,numPhysVolume1,space1,lawnum3,fullDg,2,1)
field1.stabilityParameters(beta1) # Adding stability parameters (for DG)
# // field1.setNonLocalStabilityParameters(beta1,dgnl) # Adding stability parameters (for DG)
field1.setNonLocalEqRatio(eqRatio)
# field1.gaussIntegration(0,-1,-1)
# field1.setBulkDamageBlockedMethod(1)
# field1.forceCohesiveInsertionAllIPs(bool(1),0.7)
field1.usePrimaryShapeFunction(3)
# Solver creation
# ===================================================================================
mysolver = nonLinearMechSolver(numDomain1) # Solver associated with numSolver1
# Geometry and mesh loading
# ===================================================================================
meshfile="SEN.msh" # name of mesh file
mysolver.loadModel(meshfile) # add mesh
mysolver.addDomain(field1) # add domain
# mysolver.addMaterialLaw(BulkLaw1) # add material law
# mysolver.addMaterialLaw(InterfaceLaw1)
# mysolver.addMaterialLaw(FracLaw1)
mysolver.addMaterialLaw(mylaw1)
# solver parametrisation
if Solving_mode == 1 :
mysolver.Scheme(soltype) # solver scheme
mysolver.Solver(sol) # solver choice
mysolver.snlData(nstep,ftime,tol,tolAbs) # solver parameters
mysolver.snlManageTimeStep(MaxIter,StepIncrease,StepReducFactor,NumberReduction) #timestep
# //mysolver.lineSearch(bool(0)) # lineSearch activation
# solver archiving
mysolver.stepBetweenArchiving(nstepArch) # archiving frequency for IPField
# mysolver.energyComputation(nstepArchEnergy) # archiving frequency for energy
elif Solving_mode == 2 :
mysolver.Scheme(soltype) # solver scheme
mysolver.Solver(sol) # solver choice
mysolver.snlData(nstep,ftime,tol,tolAbs) # solver parameters
mysolver.snlManageTimeStep(MaxIter,StepIncrease,StepReducFactor,NumberReduction) #timestep
mysolver.lineSearch(bool(0)) # LineSearch activation
# mysolver.explicitSpectralRadius(ftime,0.5,0.0) # Explicit integration parameters
mysolver.implicitSpectralRadius(0.05) # Implicit
# mysolver.explicitTimeStepEvaluation(nstepExpl)
# mysolver.addSystem(3,2) # adding system (numb. of var., solv. scheme)
# mysolver.addSystem(1,1)
# solver archiving
mysolver.stepBetweenArchiving(nstepArch) # archiving frequency
mysolver.energyComputation(nstepArchEnergy) # archiving frequency for energy
stiffModif = BFGSStiffnessModificationMonitoring(20, False) # stiffness will be re-computed after 20 iterations in each step
mysolver.stiffnessModification(stiffModif)
# Boundary conditions
# ===============================
# Example with constrained displacement : mysolver.displacementBC("type", numphys, comp, value)
# with: - type : Node / Edge / Face / Volume
# - numphys : number of physical entities
# - comp : number of variables on which it's applied (0 = u_x; 1 = u_y; 2 = u_z)
# - value : value at the end of the simulation
# NB: if explicit solving: displacment is replaced by a speed !
# Example with applied forces : mysolver.forceBC("type", numphys, comp, value)
if LoadingMode == 0:
#tot_disp = 0.01* 1.0*1.e-3
mysolver.displacementBC("Edge",101,0,0.) # face x = 0
mysolver.displacementBC("Edge",101,1,0.) # face x = 0
mysolver.displacementBC("Face",100,2,0.)
mysolver.forceBC("Edge",103,1, 1.0e8)
mysolver.constraintBC("Edge",103,1)
#mysolver.displacementBC("Edge",103,1,tot_disp/ftime) # face x = L
mysolver.displacementBC("Edge",103,0,0.) # face x = L
#mysolver.displacementBC("Edge",104,0,0.) # face y = 0
method = 1
mysolver.pathFollowing(True,method)
if method==0:
mysolver.setPathFollowingIncrementAdaptation(True,4,0.3)
mysolver.setPathFollowingControlType(0)
mysolver.setPathFollowingCorrectionMethod(0)
mysolver.setPathFollowingArcLengthStep(1e-1)
mysolver.setBoundsOfPathFollowingArcLengthSteps(0,0.1);
elif method==1:
mysolver.setPathFollowingIncrementAdaptationLocal(False,5,7)
mysolver.setPathFollowingLocalSteps(10.,1e-3)
mysolver.setPathFollowingLocalIncrementType(1);
mysolver.setPathFollowingSwitchCriterion(0.3)
mysolver.setBoundsOfPathFollowingLocalSteps(1e-13,4.)
elif method==2:
mysolver.setPathFollowingIncrementAdaptation(True,4,0.3)
mysolver.setPathFollowingControlType(0)
mysolver.setPathFollowingCorrectionMethod(0)
mysolver.setPathFollowingArcLengthStep(1e-1)
mysolver.setBoundsOfPathFollowingArcLengthSteps(0,2.);
# mysolver.pathFollowing(True,1)
# mysolver.setPathFollowingIncrementAdaptation(False,4)
# mysolver.setPathFollowingLocalSteps(9.0e-1,1.0e-3)
# mysolver.setPathFollowingLocalIncrementType(1);
elif LoadingMode == 1:
tot_disp = 7.0e-3 * 1.e-3
mysolver.displacementBC("Edge",101,0,0.) # face x = 0
mysolver.displacementBC("Edge",101,1,0.) # face x = 0
mysolver.displacementBC("Face",100,2,0.)
#mysolver.forceBC("Edge",103,0, 1.e8)
#mysolver.constraintBC("Edge",103,0)
mysolver.displacementBC("Edge",103,0,tot_disp/ftime) # face x = L
mysolver.displacementBC("Edge",103,1,0.) # face x = L
mysolver.displacementBC("Edge",102,1,0.)
mysolver.displacementBC("Edge",104,1,0.)
elif LoadingMode == 2:
tot_disp = 7.0e-3 * 1.e-3
mysolver.displacementBC("Edge",101,0,0.) # face x = 0
mysolver.displacementBC("Edge",101,1,0.) # face x = 0
mysolver.displacementBC("Face",100,2,0.)
mysolver.displacementBC("Edge",103,1, tot_disp/ftime)
mysolver.displacementBC("Edge",103,0, 0.)
# Variable storage
# ===============================
# Example : mysolver.internalPointBuildView("file_name", "comp", nbstep, ev)
# with: - file_name : the name of variable storage file
# - comp : saved values (IPField.SVM, IPField.SIG_XX, IPField::STRAIN_XZ,...
# see ipField.h for possibilities)
# - nbstep : storage frequency ?
# - ev : value type (default : mean)
#
# Example : mysolver.archivingForceOnPhysicalGroup("type", numphys, comp, nstep)
# with: - type : Node / Edge / Face / Volume
# - numphys : number of physical entities
# - comp : number of variables on which it's applied (0 = u_x; 1 = u_y; 2 = u_z)
# - nstep : step number between force archivingmysolver.internalPointBuildView("svm",IPField.SVM, 1, 1)
#
mysolver.internalPointBuildView("svm",IPField.SVM, 1, 1)
mysolver.internalPointBuildView("sig_xx",IPField.SIG_XX, 1, 1)
mysolver.internalPointBuildView("sig_yy",IPField.SIG_YY, 1, 1)
mysolver.internalPointBuildView("sig_zz",IPField.SIG_ZZ, 1, 1)
mysolver.internalPointBuildView("sig_xy",IPField.SIG_XY, 1, 1)
mysolver.internalPointBuildView("sig_yz",IPField.SIG_YZ, 1, 1)
mysolver.internalPointBuildView("sig_xz",IPField.SIG_XZ, 1, 1)
mysolver.internalPointBuildView("damage",IPField.DAMAGE,1,1)
# mysolver.internalPointBuildView("eff_eps",IPField.LOCAL_EQUIVALENT_STRAINS,1,1)
# mysolver.internalPointBuildView("NL_eps",IPField.NONLOCAL_EQUIVALENT_STRAINS,1,1)
# mysolver.internalPointBuildView("Blocked_element",IPField.DAMAGE_IS_BLOCKED,1,1)
mysolver.archivingForceOnPhysicalGroup("Edge", 103, 0, nstepArchForce)
mysolver.archivingForceOnPhysicalGroup("Edge", 103, 1, nstepArchForce)
mysolver.archivingForceOnPhysicalGroup("Edge", 101, 1, nstepArchForce)
mysolver.archivingForceOnPhysicalGroup("Edge", 101, 0, nstepArchForce)
mysolver.archivingNodeDisplacement(111,0, nstepArchForce)
mysolver.archivingNodeDisplacement(111,1, nstepArchForce)
# Solving
# ===========
mysolver.solve()
File deleted
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment