2#This scrip will create RVE including fibres and matrix, which can be used for all
3#three type of boundary conditions including dispalcement, traction and periodic.
6#=============================================================
7#Geometry Input parametrs
8#=============================================================
9W_wrap=0.3; H_wrap=0.1514; hgap_wrap=0.09
10W_weft=0.3; H_weft=0.0757; hgap_weft=1.2
11L_RVE=3.0; W_RVE=0.8; H_RVE=0.3;
14s_weft=hgap_weft+W_weft; s_wraf=hgap_wrap+W_wrap
15T=(H_weft/2+H_wrap/2)/2+Vgap
17L_RVE=2*s_weft; W_RVE=2*s_wraf; H_RVE=0.3;
19#=============================================================
20#Create vartices for splines to create splines
21#=============================================================
23coordx=-22.5*s_weft; coordy=T; coordz=0;
26 str1='create vertex ' + str(coordx) +' '+str(coordy)+' '+str(coordz)+' color' ; coordx=coordx+s_weft;
29 str1='create vertex ' + str(coordx) +' '+str(-coordy)+' '+str(coordz)+' color' ; coordx=coordx+s_weft;
32coordx=-22.5*s_weft; coordy=-T; coordz=s_wraf;
35 str1='create vertex ' + str(coordx) +' '+str(coordy)+' '+str(coordz)+' color' ; coordx=coordx+s_weft;
38 str1='create vertex ' + str(coordx) +' '+str(-coordy)+' '+str(coordz)+' color' ; coordx=coordx+s_weft;
41coordx=-0.5*s_weft; coordy=-T; coordz=-22*s_wraf;
44 str1='create vertex ' + str(coordx) +' '+str(coordy)+' '+str(coordz)+' color' ; coordz=coordz+s_wraf;
47 str1='create vertex ' + str(coordx) +' '+str(-coordy)+' '+str(coordz)+' color' ; coordz=coordz+s_wraf;
50coordx=0.5*s_weft; coordy=T; coordz=-22*s_wraf;
53 str1='create vertex ' + str(coordx) +' '+str(coordy)+' '+str(coordz)+' color' ; coordz=coordz+s_wraf;
56 str1='create vertex ' + str(coordx) +' '+str(-coordy)+' '+str(coordz)+' color' ; coordz=coordz+s_wraf;
59#=============================================================
60#Joint vertices to create splines
61#=============================================================
64vertices=range(1,Tvertices+1); count1=0;
66 str1='create curve spline vertex '
67 for j in range(0, 48):
68 str1=str1+' '+ str(count1+vertices[j])
69 cubit.cmd(str1); count1=count1+48;
71#=============================================================
72#Creating fibers with elliptical cross sections
73#=============================================================
75cubit.cmd('create planar surface with plane normal to curve 1 distance 0 from vertex 1')
76str1='create curve location vertex 1 direction curve 8 length '+ str(W_wrap/2); cubit.cmd(str1)
77str1='create curve location vertex 1 direction curve 7 length '+str(H_wrap/2); cubit.cmd(str1)
78cubit.cmd('create surface ellipse vertex 198 200 1')
79cubit.cmd('sweep surface 2 along curve 1')
81cubit.cmd('create planar surface with plane normal to curve 2 distance 0 from vertex 49')
82str1='create curve location vertex 49 direction curve 17 length '+ str(W_wrap/2); cubit.cmd(str1)
83str1='create curve location vertex 49 direction curve 16 length '+str(H_wrap/2); cubit.cmd(str1)
84cubit.cmd('create surface ellipse vertex 208 210 49')
85cubit.cmd('sweep surface 6 along curve 2')
87cubit.cmd('create planar surface with plane normal to curve 3 distance 0 from vertex 97')
88str1='create curve location vertex 97 direction curve 24 length '+ str(W_weft/2); cubit.cmd(str1)
89str1='create curve location vertex 97 direction curve 23 length '+str(H_weft/2); cubit.cmd(str1)
90cubit.cmd('create surface ellipse vertex 218 220 97')
91cubit.cmd('sweep surface 10 along curve 3')
93cubit.cmd('create planar surface with plane normal to curve 4 distance 0 from vertex 145')
94str1='create curve location vertex 145 direction curve 33 length '+ str(W_weft/2); cubit.cmd(str1)
95str1='create curve location vertex 145 direction curve 32 length '+str(H_weft/2); cubit.cmd(str1)
96cubit.cmd('create surface ellipse vertex 228 230 145')
97cubit.cmd('sweep surface 14 along curve 4')
99#=============================================================
100#Deletion of remaining free vertices, curves and bodies
101#=============================================================
103cubit.cmd('delete vertex all')
104cubit.cmd('delete curve all')
105cubit.cmd('delete body 7 3 5 1')
107#=============================================================
108#Create matrix and imprint and merge fibres and matrix to create common surfaces
109#=============================================================
111str1='brick x '+str(L_RVE)+' y '+str(H_RVE)+' z '+str(W_RVE); cubit.cmd(str1)
112str1='move Volume '+str(9)+' x '+str(0)+' y '+str(0)+' z '+str(s_wraf/2); cubit.cmd(str1)
114cubit.cmd('intersect volume all keep')
115cubit.cmd('delete volume 6 2 8 4')
116cubit.cmd('subtract volume 10 11 12 13 from volume 9 keep')
117cubit.cmd('delete volume 9')
118cubit.cmd('imprint volume 10 11 12 13 14')
119cubit.cmd('merge volume 10 11 12 13 14')
121#=============================================================
122#Meshing negative x, y and z face of the RVE and coping the same to its positive faces
123#=============================================================
125cubit.cmd('group 1 add surface 42 27 24')
126cubit.cmd('group "g2" add surface 44 25 28')
127cubit.cmd('group "g3" add surface 39 30 33')
128cubit.cmd('group "g4" add surface 31 40 34')
129cubit.cmd('group "g5" add surface 43')
130cubit.cmd('group "g6" add surface 41')
132cubit.cmd('surface 42 27 24 size auto factor 6')
133cubit.cmd('surface 42 27 24 scheme trimesh')
134cubit.cmd('mesh surface 42 27 24')
135cubit.cmd('surface 28 25 44 scheme mirror source surface 27 24 42 source vertex 263 target vertex 262 nosmoothing')
136cubit.cmd('mesh surface 28 25 44')
138cubit.cmd('surface 39 30 33 size auto factor 5')
139cubit.cmd('surface 39 30 33 scheme trimesh')
140cubit.cmd('mesh surface 39 30 33')
141cubit.cmd('surface 31 34 40 scheme mirror source surface 30 33 39 source vertex 257 target vertex 262 nosmoothing')
142cubit.cmd('mesh surface 31 34 40')
144cubit.cmd('surface 43 scheme trimesh')
145cubit.cmd('mesh surface 43')
146cubit.cmd('surface 41 scheme mirror source surface 43 source vertex 259 target vertex 260 nosmoothing')
147cubit.cmd('mesh surface 41')
149#=============================================================
151#=============================================================
153cubit.cmd('volume all scheme Tetmesh')
154cubit.cmd('mesh volume all')
156#=============================================================
157#Defining blocks for elastic, transversely-isotropic and potential flow problems
158#=============================================================
160vol=['14', '10,11,12,13','10', '11', '12', '13']
161mat=['MAT_ELASTIC_1','MAT_ELASTIC_TRANSISO_1','PotentialFlow_1','PotentialFlow_2','PotentialFlow_3','PotentialFlow_4']
163 cubit.cmd('set duplicate block elements on')
164 str1='block ' + str(i+1) +' volume '+vol[i]; cubit.cmd(str1)
165 str1='block ' + str(i+1) +' name "'+mat[i] + '"'; cubit.cmd(str1)
167#=============================================================
168#Material properties for matrix part
169#=============================================================
171cubit.cmd('block 1 attribute count 2')
172Em=3.5e3; Enu=0.3; #gig to mega as we used dimension in mm
174Elastic=[str(Em), str(Enu)]
176 str1='block 1 attribute index ' + str(i+1) +' '+Elastic[i]; cubit.cmd(str1)
178#=============================================================
179#Material properties for fibres
180#=============================================================
183#cubit.cmd('block 2 attribute count 5')
184#Ep=Em; Ez=Em; nup=Enu; nupz=Enu; Gzp=Em/(2*(1+Enu));
186cubit.cmd('block 2 attribute count 5')
187#Ep=40e3; Ez=230e3; nup=0.26; nupz=0.26; Gzp=24e3;
188Ep=10*Em; Ez=20*Em; nup=0.26; nupz=0.26; Gzp=5*Em;
189#Ep=2*Em; Ez=5*Em; nup=0.26; nupz=0.26; Gzp=2*Em;
191TransIso=[str(Ep), str(Ez), str(nup), str(nupz), str(Gzp)]
193 str1='block 2 attribute index ' + str(i+1) +' '+TransIso[i]; cubit.cmd(str1)
195#=============================================================
196#Material properties for interface between fibres and matrix
197#=============================================================
200cubit.cmd('set duplicate block elements on')
201str1='block 7 surface 23 26 29 32'; cubit.cmd(str1)
202str1='block 7 name "MAT_INTERF_1"'; cubit.cmd(str1)
203cubit.cmd('block 7 attribute count 4')
204str1='block 7 attribute index 1 '+str(alpha_interf); cubit.cmd(str1) #now we use 4 parameters for interface
205str1='block 7 attribute index 2 '+str(0.0); cubit.cmd(str1)
206str1='block 7 attribute index 3 '+str(0.0); cubit.cmd(str1)
207str1='block 7 attribute index 4 '+str(0.0); cubit.cmd(str1)
210#=============================================================
212#=============================================================
214Interface=['23', '26','29','32']
216 str1='sideset ' + str(i+1) +' surface '+Interface[i]; cubit.cmd(str1)
217 str1='sideset ' + str(i+1) +' name "interface'+str(i+1); cubit.cmd(str1)
219#=============================================================
220#Defining pressures for potential flow problem
221#=============================================================
223Pres=['24', '25','27','28','30','31','33','34']; count=0; count1=len(Interface);
225 str1='create pressure '+str(count+1)+' on surface '+str(Pres[count])+' magnitude 1'; cubit.cmd(str1)
226 str1='create pressure '+str(count+2)+' on surface '+str(Pres[count+1])+' magnitude -1'; cubit.cmd(str1)
227 str1='sideset '+str(count1+1)+' name "PressureIO_' + str(i+1) + '_1"'; cubit.cmd(str1)
228 str1='sideset '+str(count1+2)+' name "PressureIO_' + str(i+1) + '_2"'; cubit.cmd(str1)
229 count=count+2; count1=count1+2;
231#=============================================================
232#Defining surfaces for dispacement, traction and periodic boundary conditions
233#=============================================================
235cubit.cmd('sideset 101 surface 41 31 34 40 42 27 24') # all -ve boundary surfaces for periodic boundary conditions
236cubit.cmd('sideset 102 surface 43 30 33 39 25 44 28') # all +ve boundary surfaces for periodic boundary conditions
237cubit.cmd('sideset 103 surface 41 31 34 40 42 27 24 43 30 33 39 25 44 28') # all boundary surfaces
240#To check for indiviual side for inserstion of prisms of any problem
241#cubit.cmd('sideset 101 surface 42 27 24')
242#cubit.cmd('sideset 102 surface 25 44 28')
243#cubit.cmd('sideset 101 surface 41')
244#cubit.cmd('sideset 102 surface 43')
245#cubit.cmd('sideset 101 surface 31 34 40')
246#cubit.cmd('sideset 102 surface 30 33 39')
247#cubit.cmd('sideset 101 surface 40')
248#cubit.cmd('sideset 102 surface 39')
250#=============================================================
251#Definign zero proessrues for potential flow problem (This should be of the same order as PotentialFlow blocks )
252#=============================================================
254zeroPressureNode=[140, 106, 279, 272]
256 str1='nodeset ' + str(i+1) + ' node ' + str(zeroPressureNode[i]); cubit.cmd(str1)
257 str1='nodeset ' + str(i+1)+' name "ZeroPressure_' + str(i+1)+ '"'; cubit.cmd(str1)
259#=============================================================
260#Saving input RVE file
261#=============================================================
263cubit.cmd('save as "/Users/zahur/Documents/moFEM/mofem-cephas/mofem_v0.2/users_modules/homogenisation/meshes/Textile_RVE_TransIso.cub" overwrite')