v0.14.0
rve_mech_plas_interface_periodic.cpp
Go to the documentation of this file.
1 /** \file small_strain_plasticity.cpp
2  * \ingroup small_strain_plasticity
3  * \brief Small strain plasticity example
4  *
5  */
6 
7 /* This file is part of MoFEM.
8  * MoFEM is free software: you can redistribute it and/or modify it under
9  * the terms of the GNU Lesser General Public License as published by the
10  * Free Software Foundation, either version 3 of the License, or (at your
11  * option) any later version.
12  *
13  * MoFEM is distributed in the hope that it will be useful, but WITHOUT
14  * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
15  * FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
16  * License for more details.
17  *
18  * You should have received a copy of the GNU Lesser General Public
19  * License along with MoFEM. If not, see <http://www.gnu.org/licenses/>. */
20 
21 
22 #include <MoFEM.hpp>
23 using namespace MoFEM;
24 
25 #include <boost/program_options.hpp>
26 using namespace std;
27 namespace po = boost::program_options;
28 
29 #include <boost/numeric/ublas/vector_proxy.hpp>
30 #include <boost/numeric/ublas/matrix.hpp>
31 #include <boost/numeric/ublas/matrix_proxy.hpp>
32 #include <boost/numeric/ublas/vector.hpp>
33 #include <boost/numeric/ublas/symmetric.hpp>
34 
35 
36 #include <MethodForForceScaling.hpp>
37 #include <TimeForceScale.hpp>
38 #include <VolumeCalculation.hpp>
39 
40 #include <VolumeCalculation.hpp>
41 #include <BCs_RVELagrange_Disp.hpp>
42 #include <BCs_RVELagrange_Trac.hpp>
43 #include "BCs_RVELagrange_Periodic.hpp"
44 #include <BCs_RVELagrange_Trac_Rigid_Trans.hpp>
45 
46 #include <CohesiveInterfaceElement.hpp>
47 using namespace CohesiveElement;
48 
49 #ifndef WITH_ADOL_C
50 #error "MoFEM need to be compiled with ADOL-C"
51 #endif
52 #include <adolc/adolc.h>
53 #include <NonLinearElasticElement.hpp>
54 #include <SmallTransverselyIsotropic.hpp>
55 
56 #include <adolc/adolc.h>
57 #include <SmallStrainPlasticity.hpp>
58 #include <SmallStrainPlasticityMaterialModels.hpp>
59 
60 
61 #include <Projection10NodeCoordsOnField.hpp>
62 #include <PostProcOnRefMesh.hpp>
63 
64 #include <string>
65 
66 using namespace boost::numeric;
67 
68 
69 
70 
71 static char help[] =
72 "...\n"
73 "\n";
74 
75 
76 //=================================================================================================================================
77 //Define class and multindex container to store data for traiangles on the boundary of the RVE (it cannot be defined within main)
78 //=================================================================================================================================
79 
80 struct Face_CenPos_Handle {
81  double xcoord, ycoord, zcoord;
82  const EntityHandle Tri_Hand;
83  Face_CenPos_Handle(double _xcoord, double _ycoord, double _zcoord, const EntityHandle _Tri_Hand):xcoord(_xcoord),
84  ycoord(_ycoord), zcoord(_zcoord), Tri_Hand(_Tri_Hand) {}
85 };
86 
87 struct xcoord_tag {}; //tags to used in the multindex container
88 struct ycoord_tag {};
89 struct zcoord_tag {};
90 struct Tri_Hand_tag {};
91 struct Composite_xyzcoord {};
92 
93 typedef multi_index_container<
94 Face_CenPos_Handle,
95 indexed_by<
96 ordered_non_unique<
97 tag<xcoord_tag>, member<Face_CenPos_Handle,double,&Face_CenPos_Handle::xcoord> >,
98 
99 ordered_non_unique<
100 tag<ycoord_tag>, member<Face_CenPos_Handle,double,&Face_CenPos_Handle::ycoord> >,
101 
102 ordered_non_unique<
103 tag<zcoord_tag>, member<Face_CenPos_Handle,double,&Face_CenPos_Handle::zcoord> >,
104 
105 ordered_unique<
106 tag<Tri_Hand_tag>, member<Face_CenPos_Handle,const EntityHandle,&Face_CenPos_Handle::Tri_Hand> >,
107 
108 ordered_unique<
109 tag<Composite_xyzcoord>,
110 composite_key<
111 Face_CenPos_Handle,
112 member<Face_CenPos_Handle,double,&Face_CenPos_Handle::xcoord>,
113 member<Face_CenPos_Handle,double,&Face_CenPos_Handle::ycoord>,
114 member<Face_CenPos_Handle,double,&Face_CenPos_Handle::zcoord> > >
115 > > Face_CenPos_Handle_multiIndex;
116 
117 //============================================================================================
118 //============================================================================================
119 
120 
121 
122 int main(int argc, char *argv[]) {
123  MoFEM::Core::Initialize(&argc,&argv,(char *)0,help);
124 
125  try {
126 
127  moab::Core mb_instance;
128  moab::Interface& moab = mb_instance;
129  ParallelComm* pcomm = ParallelComm::get_pcomm(&moab,MYPCOMM_INDEX);
130  if(pcomm == NULL) pcomm = new ParallelComm(&moab,PETSC_COMM_WORLD);
131 
132  PetscBool flg = PETSC_TRUE;
133  char mesh_file_name[255];
134  ierr = PetscOptionsGetString(PETSC_NULL,"-my_file",mesh_file_name,255,&flg); CHKERRQ(ierr);
135  if(flg != PETSC_TRUE) {
136  SETERRQ(PETSC_COMM_SELF,1,"*** ERROR -my_file (MESH FILE NEEDED)");
137  }
138 
139  PetscInt order;
140  ierr = PetscOptionsGetInt(PETSC_NULL,"-my_order",&order,&flg); CHKERRQ(ierr);
141  if(flg != PETSC_TRUE) {
142  order = 2;
143  }
144  PetscInt bubble_order;
145  ierr = PetscOptionsGetInt(PETSC_NULL,"-my_bubble_order",&bubble_order,&flg); CHKERRQ(ierr);
146  if(flg != PETSC_TRUE) {
147  bubble_order = order;
148  }
149 
150  // use this if your mesh is partitioned and you run code on parts,
151  // you can solve very big problems
152  PetscBool is_partitioned = PETSC_FALSE;
153  ierr = PetscOptionsGetBool(PETSC_NULL,"-my_is_partitioned",&is_partitioned,&flg); CHKERRQ(ierr);
154 
155  //Applied strain on the RVE (vector of length 6) strain=[xx, yy, zz, xy, xz, zy]^T
156  double mygiven_strain[6];
157  int nmax=6;
158  ierr = PetscOptionsGetRealArray(PETSC_NULL,"-my_given_strain",mygiven_strain,&nmax,&flg); CHKERRQ(ierr);
159  VectorDouble given_strain;
160  given_strain.resize(6);
161  cblas_dcopy(6, &mygiven_strain[0], 1, &given_strain(0), 1);
162  cout<<"given_strain ="<<given_strain<<endl;
163 
164  //Read mesh to MOAB
165 
166  if(is_partitioned == PETSC_TRUE) {
167  const char *option;
168  option = "PARALLEL=BCAST_DELETE;PARALLEL_RESOLVE_SHARED_ENTS;PARTITION=PARALLEL_PARTITION;";
169  rval = moab.load_file(mesh_file_name, 0, option); CHKERRQ_MOAB(rval);
170  } else {
171  const char *option;
172  option = "";
173  rval = moab.load_file(mesh_file_name, 0, option); CHKERRQ_MOAB(rval);
174  }
175 
176 
177 
178  double young_modulus = 1;
179  SmallStrainParaboloidalPlasticity cp;
180  {
181  cp.tAgs.resize(3);
182  cp.tAgs[0] = 3;
183  cp.tAgs[1] = 4;
184  cp.tAgs[2] = 5;
185  cp.tOl = 1e-12;
186 
187  double poisson_ratio = 0.25;
188  cp.sIgma_yt = 1;
189  cp.sIgma_yc = 1;
190 
191  cp.Ht = 0.1;
192  cp.Hc = 0.1;
193 
194  cp.nup = 0.3;
195 
196 
197  {
198  ierr = PetscOptionsGetReal(0,"-my_young_modulus",&young_modulus,0); CHKERRQ(ierr);
199  ierr = PetscOptionsGetReal(0,"-my_poisson_ratio",&poisson_ratio,0); CHKERRQ(ierr);
200  cp.mu = MU(young_modulus,poisson_ratio);
201  cp.lambda = LAMBDA(young_modulus,poisson_ratio);
202  ierr = PetscOptionsGetReal(0,"-my_sigma_yt",&cp.sIgma_yt,0); CHKERRQ(ierr);
203  ierr = PetscOptionsGetReal(0,"-my_sigma_yc",&cp.sIgma_yc,0); CHKERRQ(ierr);
204 
205  ierr = PetscOptionsGetReal(0,"-my_Ht",&cp.Ht,0); CHKERRQ(ierr);
206  ierr = PetscOptionsGetReal(0,"-my_Hc",&cp.Hc,0); CHKERRQ(ierr);
207 
208  ierr = PetscOptionsGetReal(0,"-my_nt",&cp.nt,0); CHKERRQ(ierr);
209  ierr = PetscOptionsGetReal(0,"-my_nc",&cp.nc,0); CHKERRQ(ierr);
210 
211 
212  ierr = PetscOptionsGetReal(0,"-my_nup",&cp.nup,0); CHKERRQ(ierr);
213  }
214 
215  PetscPrintf(PETSC_COMM_WORLD,"young_modulus = %4.2e \n",young_modulus);
216  PetscPrintf(PETSC_COMM_WORLD,"poisson_ratio = %4.2e \n",poisson_ratio);
217 
218  PetscPrintf(PETSC_COMM_WORLD,"sIgma_yt = %4.2e \n",cp.sIgma_yt);
219  PetscPrintf(PETSC_COMM_WORLD,"sIgma_yc = %4.2e \n",cp.sIgma_yc);
220 
221  PetscPrintf(PETSC_COMM_WORLD,"nup = %4.2e \n",cp.nup);
222 
223  cp.sTrain.resize(6,false);
224  cp.sTrain.clear();
225  cp.plasticStrain.resize(6,false);
226  cp.plasticStrain.clear();
227  cp.internalVariables.resize(2,false);
228  cp.internalVariables.clear();
229  cp.createMatAVecR();
230  cp.snesCreate();
231  // ierr = SNESSetFromOptions(cp.sNes); CHKERRQ(ierr);
232  }
233 
234 
235  // double young_modulus = 1;
236  // SmallStrainJ2Plasticity cp;
237  // {
238  // cp.tAgs.resize(3);
239  // cp.tAgs[0] = 3;
240  // cp.tAgs[1] = 4;
241  // cp.tAgs[2] = 5;
242  // cp.tOl = 1e-12;
243  //
244  // double poisson_ratio = 0.25;
245  // cp.sIgma_y = 1;
246  // cp.H = 0.1;
247  // cp.K = 0;
248  // cp.phi = 1;
249  // {
250  // ierr = PetscOptionsGetReal(0,"-my_young_modulus",&young_modulus,0); CHKERRQ(ierr);
251  // ierr = PetscOptionsGetReal(0,"-my_poisson_ratio",&poisson_ratio,0); CHKERRQ(ierr);
252  // cp.mu = MU(young_modulus,poisson_ratio);
253  // cp.lambda = LAMBDA(young_modulus,poisson_ratio);
254  // ierr = PetscOptionsGetReal(0,"-my_sigma_y",&cp.sIgma_y,0); CHKERRQ(ierr);
255  // ierr = PetscOptionsGetReal(0,"-my_H",&cp.H,0); CHKERRQ(ierr);
256  // ierr = PetscOptionsGetReal(0,"-my_K",&cp.K,0); CHKERRQ(ierr);
257  // ierr = PetscOptionsGetReal(0,"-my_phi",&cp.phi,0); CHKERRQ(ierr);
258  // }
259  //
260  // PetscPrintf(PETSC_COMM_WORLD,"young_modulus = %4.2e \n",young_modulus);
261  // PetscPrintf(PETSC_COMM_WORLD,"poisson_ratio = %4.2e \n",poisson_ratio);
262  // PetscPrintf(PETSC_COMM_WORLD,"sIgma_y = %4.2e \n",cp.sIgma_y);
263  // PetscPrintf(PETSC_COMM_WORLD,"H = %4.2e \n",cp.H);
264  // PetscPrintf(PETSC_COMM_WORLD,"K = %4.2e \n",cp.K);
265  // PetscPrintf(PETSC_COMM_WORLD,"phi = %4.2e \n",cp.phi);
266  //
267  //
268  // cp.sTrain.resize(6,false);
269  // cp.sTrain.clear();
270  // cp.plasticStrain.resize(6,false);
271  // cp.plasticStrain.clear();
272  // cp.internalVariables.resize(7,false);
273  // cp.internalVariables.clear();
274  // cp.createMatAVecR();
275  // cp.snesCreate();
276  // // ierr = SNESSetFromOptions(cp.sNes); CHKERRQ(ierr);
277  // }
278 
279 
280  //Create MoFEM (Joseph) database
281  MoFEM::Core core(moab);
282  MoFEM::Interface& m_field = core;
283 
284  vector<BitRefLevel> bit_levels;
285  Tag th_meshset_info;
286  int def_meshset_info[2] = {0,0};
287  rval = moab.tag_get_handle("MESHSET_INFO",2,MB_TYPE_INTEGER,th_meshset_info,MB_TAG_CREAT|MB_TAG_SPARSE,&def_meshset_info);
288  int meshset_data[2];
289  EntityHandle root = moab.get_root_set();
290  rval = moab.tag_get_data(th_meshset_info,&root,1,meshset_data); CHKERRQ_MOAB(rval);
291  if(meshset_data[0]==0) {
292  meshset_data[0] = 1;
293  rval = moab.tag_set_data(th_meshset_info,&root,1,meshset_data); CHKERRQ_MOAB(rval);
294  }
295  bit_levels.push_back(BitRefLevel().set(meshset_data[0]-1));
296 
297 
298 // EntityHandle out_meshset;
299 // rval = moab.create_meshset(MESHSET_SET,out_meshset); CHKERRQ_MOAB(rval);
300 // // ierr = m_field.get_problem_finite_elements_entities("POTENTIAL_PROBLEM","POTENTIAL_ELEM",out_meshset); CHKERRQ(ierr);
301 // ierr = m_field.get_entities_by_ref_level(bit_levels.back(),BitRefLevel().set(),out_meshset); CHKERRQ(ierr);
302 // Range LatestRefinedTets;
303 // rval = moab.get_entities_by_type(out_meshset, MBTET,LatestRefinedTets,true); CHKERRQ_MOAB(rval);
304 // Range LatestRefinedPrisms;
305 // rval = moab.get_entities_by_type(out_meshset, MBPRISM,LatestRefinedPrisms,true); CHKERRQ_MOAB(rval);
306 
307 
308  BitRefLevel problem_bit_level = bit_levels.back();
309 
310  // const clock_t begin_time = clock();
311  ierr = m_field.build_fields(); CHKERRQ(ierr);
312  ierr = m_field.build_finite_elements(); CHKERRQ(ierr);
313 
314 
315  Range interface_prims_on_problem_bit_level;
316  ierr = m_field.get_entities_by_type_and_ref_level(
317  problem_bit_level,BitRefLevel().set(),MBPRISM,interface_prims_on_problem_bit_level
318  ); CHKERRQ(ierr);
319  Range tets_on_problem_bit_level;
320  ierr = m_field.get_entities_by_type_and_ref_level(
321  problem_bit_level,BitRefLevel().set(),MBTET,tets_on_problem_bit_level
322  ); CHKERRQ(ierr);
323 
324  cout<<"interface_prims_on_problem_bit_level "<<interface_prims_on_problem_bit_level.size()<<endl;
325  cout<<"tets_on_problem_bit_level "<<tets_on_problem_bit_level.size()<<endl;
326 
327  //=======================================================================================================
328  //Add Periodic Prisims Between Triangles on -ve and +ve faces to implement periodic bounary conditions
329  //=======================================================================================================
330  Range preriodic_prisms;
331  //Populating the Multi-index container with -ve triangles
332 
333 
334  Range SurTrisNeg; //Old + new traiangles on side set 101
335  ierr = m_field.get_cubit_msId_entities_by_dimension(101,SIDESET,2,SurTrisNeg,true); CHKERRQ(ierr);
336  ierr = PetscPrintf(PETSC_COMM_WORLD,"Old+New triangles on SideSet 101 = %d\n",SurTrisNeg.size()); CHKERRQ(ierr);
337 
338  Range AllTris; //All new traiangles in the new-bitlevel
339  ierr = m_field.get_entities_by_type_and_ref_level(problem_bit_level,BitRefLevel().set(),MBTRI,AllTris); CHKERRQ(ierr);
340  ierr = PetscPrintf(PETSC_COMM_WORLD,"All triangles in the new bit-level = %d\n",AllTris.size()); CHKERRQ(ierr);
341 
342  //Only triangles in the new bit-level
343  SurTrisNeg = intersect(AllTris,SurTrisNeg);
344  ierr = PetscPrintf(PETSC_COMM_WORLD,"New triangles on SideSet 101 = %d\n",SurTrisNeg.size()); CHKERRQ(ierr);
345 
346 // cout<<"SurTrisNeg "<<SurTrisNeg<<endl;
347  Face_CenPos_Handle_multiIndex Face_CenPos_Handle_varNeg, Face_CenPos_Handle_varPos;
348  double TriCen[3], coords_Tri[9];
349 
350 
351 
352  double roundfact=1e3;
353 
354  for(Range::iterator it = SurTrisNeg.begin(); it!=SurTrisNeg.end(); it++) {
355  const EntityHandle* conn_face; int num_nodes_Tri;
356 
357  //get nodes attached to the face
358  rval = moab.get_connectivity(*it,conn_face,num_nodes_Tri,true); CHKERRQ_MOAB(rval);
359  //get nodal coordinates
360  rval = moab.get_coords(conn_face,num_nodes_Tri,coords_Tri); CHKERRQ_MOAB(rval);
361 
362  //Find triangle centriod
363  TriCen[0]= (coords_Tri[0]+coords_Tri[3]+coords_Tri[6])/3.0;
364  TriCen[1]= (coords_Tri[1]+coords_Tri[4]+coords_Tri[7])/3.0;
365  TriCen[2]= (coords_Tri[2]+coords_Tri[5]+coords_Tri[8])/3.0;
366 
367 // cout<<"\nTriCen[0]= "<<TriCen[0] << " TriCen[1]= "<< TriCen[1] << " TriCen[2]= "<< TriCen[2] <<endl;
368  TriCen[0]=round(TriCen[0]*roundfact)/roundfact;
369  TriCen[1]=round(TriCen[1]*roundfact)/roundfact;
370  TriCen[2]=round(TriCen[2]*roundfact)/roundfact;
371 
372 // cout<<"TriCen[0]= "<<TriCen[0] << " TriCen[1]= "<< TriCen[1] << " TriCen[2]= "<< TriCen[2] <<endl;
373  //fill the multi-index container with centriod coordinates and triangle handles
374  Face_CenPos_Handle_varNeg.insert(Face_CenPos_Handle(TriCen[0], TriCen[1], TriCen[2], *it));
375  // for(int ii=0; ii<3; ii++) cout<<"TriCen "<<TriCen[ii]<<endl;
376  // cout<<endl<<endl;
377  }
378 
379  // double aaa;
380  // aaa=0.5011;
381  // cout<<"\n\n\n\n\nfloor(aaa+0.5) = "<<double(int(aaa*roundfact+0.5))/roundfact<<endl<<endl<<endl<<endl;
382  // aaa=-0.5011;
383  // cout<<"\n\n\n\n\nfloor(aaa+0.5) = "<<double(int(aaa*roundfact-0.5))/roundfact<<endl<<endl<<endl<<endl;
384 
385 
386  //Populating the Multi-index container with +ve triangles
387  Range SurTrisPos;
388  ierr = m_field.get_cubit_msId_entities_by_dimension(102,SIDESET,2,SurTrisPos,true); CHKERRQ(ierr);
389  ierr = PetscPrintf(PETSC_COMM_WORLD,"Old+New triangles on SideSet 102 = %d\n",SurTrisPos.size()); CHKERRQ(ierr);
390 
391  //Only triangles in the new bit-level
392  SurTrisPos = intersect(AllTris,SurTrisPos);
393  ierr = PetscPrintf(PETSC_COMM_WORLD,"New triangles on SideSet 102 = %d\n",SurTrisPos.size()); CHKERRQ(ierr);
394 
395 
396 
397 
398  for(Range::iterator it = SurTrisPos.begin(); it!=SurTrisPos.end(); it++) {
399  const EntityHandle* conn_face; int num_nodes_Tri;
400 
401  //get nodes attached to the face
402  rval = moab.get_connectivity(*it,conn_face,num_nodes_Tri,true); CHKERRQ_MOAB(rval);
403  //get nodal coordinates
404  rval = moab.get_coords(conn_face,num_nodes_Tri,coords_Tri); CHKERRQ_MOAB(rval);
405 
406  //Find triangle centriod
407  TriCen[0]= (coords_Tri[0]+coords_Tri[3]+coords_Tri[6])/3.0;
408  TriCen[1]= (coords_Tri[1]+coords_Tri[4]+coords_Tri[7])/3.0;
409  TriCen[2]= (coords_Tri[2]+coords_Tri[5]+coords_Tri[8])/3.0;
410 
411  //round values to 3 disimal places
412  TriCen[0]=round(TriCen[0]*roundfact)/roundfact;
413  TriCen[1]=round(TriCen[1]*roundfact)/roundfact;
414  TriCen[2]=round(TriCen[2]*roundfact)/roundfact;
415 
416 // cout<<"TriCen[0]= "<<TriCen[0] << " TriCen[1]= "<< TriCen[1] << " TriCen[2]= "<< TriCen[2] <<endl;
417 
418  //fill the multi-index container with centriod coordinates and triangle handles
419  Face_CenPos_Handle_varPos.insert(Face_CenPos_Handle(TriCen[0], TriCen[1], TriCen[2], *it));
420  }
421 
422  //Find minimum and maximum X, Y and Z coordinates of the RVE (using multi-index container)
423  double XcoordMin, YcoordMin, ZcoordMin, XcoordMax, YcoordMax, ZcoordMax;
424  typedef Face_CenPos_Handle_multiIndex::index<xcoord_tag>::type::iterator Tri_Xcoord_iterator;
425  typedef Face_CenPos_Handle_multiIndex::index<ycoord_tag>::type::iterator Tri_Ycoord_iterator;
426  typedef Face_CenPos_Handle_multiIndex::index<zcoord_tag>::type::iterator Tri_Zcoord_iterator;
427  Tri_Xcoord_iterator XcoordMin_it, XcoordMax_it;
428  Tri_Ycoord_iterator YcoordMin_it, YcoordMax_it;
429  Tri_Zcoord_iterator ZcoordMin_it, ZcoordMax_it;
430 
431  //XcoordMax_it-- because .end() will point iterator after the data range but .begin() will point the iteratore to the first value of range
432  XcoordMin_it=Face_CenPos_Handle_varNeg.get<xcoord_tag>().begin(); XcoordMin=XcoordMin_it->xcoord;
433  XcoordMax_it=Face_CenPos_Handle_varPos.get<xcoord_tag>().end(); XcoordMax_it--; XcoordMax=XcoordMax_it->xcoord;
434  YcoordMin_it=Face_CenPos_Handle_varNeg.get<ycoord_tag>().begin(); YcoordMin=YcoordMin_it->ycoord;
435  YcoordMax_it=Face_CenPos_Handle_varPos.get<ycoord_tag>().end(); YcoordMax_it--; YcoordMax=YcoordMax_it->ycoord;
436  ZcoordMin_it=Face_CenPos_Handle_varNeg.get<zcoord_tag>().begin(); ZcoordMin=ZcoordMin_it->zcoord;
437  ZcoordMax_it=Face_CenPos_Handle_varPos.get<zcoord_tag>().end(); ZcoordMax_it--; ZcoordMax=ZcoordMax_it->zcoord;
438 
439  cout<<"XcoordMin "<<XcoordMin << " XcoordMax "<<XcoordMax <<endl;
440  cout<<"YcoordMin "<<YcoordMin << " YcoordMax "<<YcoordMax <<endl;
441  cout<<"ZcoordMin "<<ZcoordMin << " ZcoordMax "<<ZcoordMax <<endl;
442 
443  /*double LxRVE, LyRVE, LzRVE;
444  LxRVE=XcoordMax-XcoordMin;
445  LyRVE=YcoordMax-YcoordMin;
446  LzRVE=ZcoordMax-ZcoordMin;*/
447 
448  //Creating Prisims between triangles on -ve and +ve faces
449  typedef Face_CenPos_Handle_multiIndex::index<Tri_Hand_tag>::type::iterator Tri_Hand_iterator;
450  Tri_Hand_iterator Tri_Neg;
451  typedef Face_CenPos_Handle_multiIndex::index<Composite_xyzcoord>::type::iterator xyzcoord_iterator;
452  xyzcoord_iterator Tri_Pos;
453  double XPos, YPos, ZPos;
454  //int count=0;
455 
456  //loop over -ve triangles to create prisims elemenet between +ve and -ve triangles
457  int count=1;
458  int count1=0;
459  for(Range::iterator it = SurTrisNeg.begin(); it!=SurTrisNeg.end(); it++) {
460 
461 // count1++;
462 // cout<<"count1 = "<<count1<< endl;
463 // if(count1 == 24) {
464 
465 
466  Tri_Neg=Face_CenPos_Handle_varNeg.get<Tri_Hand_tag>().find(*it);
467  //cout<<"xcoord= "<<Tri_iit->xcoord << " ycoord= "<< Tri_iit->ycoord << " ycoord= "<< Tri_iit->zcoord <<endl;
468 
469  //corresponding +ve triangle
470  if(Tri_Neg->xcoord==XcoordMin){XPos=XcoordMax; YPos=Tri_Neg->ycoord; ZPos=Tri_Neg->zcoord;};
471  if(Tri_Neg->ycoord==YcoordMin){XPos=YPos=Tri_Neg->xcoord; YPos=YcoordMax; ZPos=Tri_Neg->zcoord;};
472  if(Tri_Neg->zcoord==ZcoordMin){XPos=YPos=Tri_Neg->xcoord; YPos=Tri_Neg->ycoord; ZPos=ZcoordMax; };
473  Tri_Pos=Face_CenPos_Handle_varPos.get<Composite_xyzcoord>().find(boost::make_tuple(XPos, YPos, ZPos));
474 
475 // cout<<"Tri_Neg->xcoord= "<<Tri_Neg->xcoord << " Tri_Neg->ycoord "<< Tri_Neg->ycoord << " Tri_Neg->zcoord= "<< Tri_Neg->zcoord <<endl;
476 // cout<<"Tri_Pos->xcoord= "<<Tri_Pos->xcoord << " Tri_Pos->ycoord "<< Tri_Pos->ycoord << " Tri_Pos->zcoord= "<< Tri_Pos->zcoord <<endl;
477 
478 
479  //+ve and -ve nodes and their coords (+ve and -ve tiangles nodes can have matching problems, which can produce twisted prism)
480  EntityHandle PrismNodes[6];
481  vector<EntityHandle> TriNodesNeg, TriNodesPos;
482  double CoordNodeNeg[9], CoordNodePos[9];
483  rval = moab.get_connectivity(&(Tri_Neg->Tri_Hand),1,TriNodesNeg,true); CHKERRQ_MOAB(rval);
484  rval = moab.get_connectivity(&(Tri_Pos->Tri_Hand),1,TriNodesPos,true); CHKERRQ_MOAB(rval);
485  rval = moab.get_coords(&TriNodesNeg[0],3,CoordNodeNeg); MOAB_THROW(rval);
486  rval = moab.get_coords(&TriNodesPos[0],3,CoordNodePos); MOAB_THROW(rval);
487  for(int ii=0; ii<3; ii++){
488  PrismNodes[ii]=TriNodesNeg[ii];
489  }
490 
491 // for(int ii=0; ii<3; ii++){
492 // cout<<"xcoord= "<<CoordNodeNeg[3*ii] << " ycoord= "<< CoordNodeNeg[3*ii+1] << " zcoord= "<< CoordNodeNeg[3*ii+2] <<endl;
493 // }
494 // for(int ii=0; ii<3; ii++){
495 // cout<<"xcoord= "<<CoordNodePos[3*ii] << " ycoord= "<< CoordNodePos[3*ii+1] << " zcoord= "<< CoordNodePos[3*ii+2] <<endl;
496 // }
497 
498  //Match exact nodes to each other to avoide the problem of twisted prisms
499  double XNodeNeg, YNodeNeg, ZNodeNeg, XNodePos, YNodePos, ZNodePos;
500  for(int ii=0; ii<3; ii++){
501  if(Tri_Neg->xcoord==XcoordMin){XNodeNeg=XcoordMax; YNodeNeg=CoordNodeNeg[3*ii+1]; ZNodeNeg=CoordNodeNeg[3*ii+2];};
502  if(Tri_Neg->ycoord==YcoordMin){XNodeNeg=CoordNodeNeg[3*ii]; YNodeNeg=YcoordMax; ZNodeNeg=CoordNodeNeg[3*ii+2];};
503  if(Tri_Neg->zcoord==ZcoordMin){XNodeNeg=CoordNodeNeg[3*ii]; YNodeNeg=CoordNodeNeg[3*ii+1]; ZNodeNeg=ZcoordMax;};
504  for(int jj=0; jj<3; jj++){
505  //Round nodal coordinates to 3 dicimal places only for comparison
506  XNodeNeg=round(XNodeNeg*roundfact)/roundfact;
507  YNodeNeg=round(YNodeNeg*roundfact)/roundfact;
508  ZNodeNeg=round(ZNodeNeg*roundfact)/roundfact;
509 
510  XNodePos=CoordNodePos[3*jj]; YNodePos=CoordNodePos[3*jj+1]; ZNodePos=CoordNodePos[3*jj+2];
511  XNodePos=round(XNodePos*roundfact)/roundfact;
512  YNodePos=round(YNodePos*roundfact)/roundfact;
513  ZNodePos=round(ZNodePos*roundfact)/roundfact;
514 
515  if(XNodeNeg==XNodePos && YNodeNeg==YNodePos && ZNodeNeg==ZNodePos){
516  PrismNodes[3+ii]=TriNodesPos[jj];
517  break;
518  }
519  }
520  }
521  //prism nodes and their coordinates
522  double CoordNodesPrisms[18];
523  rval = moab.get_coords(PrismNodes,6,CoordNodesPrisms); MOAB_THROW(rval);
524 // for(int ii=0; ii<6; ii++){
525 // cout<<"xcoord= "<<CoordNodesPrisms[3*ii] << " ycoord= "<< CoordNodesPrisms[3*ii+1] << " zcoord= "<< CoordNodesPrisms[3*ii+2] <<endl;
526 // }
527 // cout<<endl<<endl;
528 
529  //insertion of individula prism element and its addition to range preriodic_prisms
530  EntityHandle PeriodicPrism;
531  rval = moab.create_element(MBPRISM,PrismNodes,6,PeriodicPrism); CHKERRQ_MOAB(rval);
532  preriodic_prisms.insert(PeriodicPrism);
533 
534 // //to see individual prisms
535 // Range Prism1;
536 // Prism1.insert(PeriodicPrism);
537 // EntityHandle out_meshset1;
538 // rval = moab.create_meshset(MESHSET_SET,out_meshset1); CHKERRQ_MOAB(rval);
539 // rval = moab.add_entities(out_meshset1,Prism1); CHKERRQ_MOAB(rval);
540 // ostringstream sss;
541 // sss << "Prism" << count << ".vtk"; count++;
542 // rval = moab.write_file(sss.str().c_str(),"VTK","",&out_meshset1,1); CHKERRQ_MOAB(rval);
543 //
544 // string aaa;
545 // cin>>aaa;
546 //
547 // }
548 
549  }
550 
551 
552 
553 
554  //Adding the newely added prisms to the bit-level (which is known as seeding)
555  /*{
556  Range ents;
557  moab.get_adjacencies(preriodic_prisms,1,true,ents,moab::Interface::UNION);
558  moab.get_adjacencies(preriodic_prisms,2,true,ents,moab::Interface::UNION);
559  ierr = m_field.seed_ref_level(ents.subset_by_type(MBEDGE),problem_bit_level); CHKERRQ(ierr);
560  ierr = m_field.seed_ref_level(ents.subset_by_type(MBQUAD),problem_bit_level); CHKERRQ(ierr);
561  }*/
562  ierr = m_field.seed_ref_level(preriodic_prisms,problem_bit_level); CHKERRQ(ierr);
563  ierr = m_field.seed_finite_elements(preriodic_prisms); CHKERRQ(ierr);
564 
565 
566  //=======================================================================================================
567  //=======================================================================================================
568  //=======================================================================================================
569 
570  EntityHandle out_meshset;
571  rval = moab.create_meshset(MESHSET_SET,out_meshset); CHKERRQ_MOAB(rval);
572  ierr = m_field.get_entities_by_ref_level(bit_levels.back(),BitRefLevel().set(),out_meshset); CHKERRQ(ierr);
573 
574  Range LatestRefinedTets;
575  rval = moab.get_entities_by_type(out_meshset, MBTET,LatestRefinedTets,true); CHKERRQ_MOAB(rval);
576 
577  Range Interface_Periodic_Prisms;
578  rval = moab.get_entities_by_type(out_meshset, MBPRISM,Interface_Periodic_Prisms,true); CHKERRQ_MOAB(rval);
579  cout<<"Interface_Periodic_Prisms " <<Interface_Periodic_Prisms.size()<<endl;
580 
581 
582  //Fields
583  ierr = m_field.add_field("DISPLACEMENT",H1,AINSWORTH_LEGENDRE_BASE,3); CHKERRQ(ierr);
584  ierr = m_field.add_field("MESH_NODE_POSITIONS",H1,AINSWORTH_LEGENDRE_BASE,3,MB_TAG_SPARSE,MF_ZERO); CHKERRQ(ierr);
585  ierr = m_field.add_field("LAGRANGE_MUL_PERIODIC",H1,AINSWORTH_LEGENDRE_BASE,3); CHKERRQ(ierr); //For lagrange multipliers to control the periodic motion
586  ierr = m_field.add_field("LAGRANGE_MUL_RIGID_TRANS",NOFIELD,NOBASE,3); CHKERRQ(ierr); //To control the rigid body motion (3 Traslations and 3
587 
588  //add entitities (by tets) to the field
589  ierr = m_field.add_ents_to_field_by_type(0,MBTET,"DISPLACEMENT"); CHKERRQ(ierr);
590  ierr = m_field.add_ents_to_field_by_type(0,MBTET,"MESH_NODE_POSITIONS"); CHKERRQ(ierr);
591 
592 // ierr = PetscPrintf(PETSC_COMM_WORLD,"number of triangle in SideSet 101 = %d\n",SurTrisNeg.size()); CHKERRQ(ierr);
593  ierr = m_field.add_ents_to_field_by_type(SurTrisNeg,MBTRI,"LAGRANGE_MUL_PERIODIC"); CHKERRQ(ierr);
594 
595  //set app. order
596  ierr = m_field.set_field_order(0,MBTET,"DISPLACEMENT",bubble_order); CHKERRQ(ierr);
597  ierr = m_field.set_field_order(0,MBTRI,"DISPLACEMENT",order); CHKERRQ(ierr);
598  ierr = m_field.set_field_order(0,MBEDGE,"DISPLACEMENT",order); CHKERRQ(ierr);
599  ierr = m_field.set_field_order(0,MBVERTEX,"DISPLACEMENT",1); CHKERRQ(ierr);
600 
601  ierr = m_field.set_field_order(0,MBTET,"MESH_NODE_POSITIONS",order>1 ? 2 : 1); CHKERRQ(ierr);
602  ierr = m_field.set_field_order(0,MBTRI,"MESH_NODE_POSITIONS",order>1 ? 2 : 1); CHKERRQ(ierr);
603  ierr = m_field.set_field_order(0,MBEDGE,"MESH_NODE_POSITIONS",order>1 ? 2 : 1); CHKERRQ(ierr);
604  ierr = m_field.set_field_order(0,MBVERTEX,"MESH_NODE_POSITIONS",1); CHKERRQ(ierr);
605 
606  ierr = m_field.set_field_order(0,MBTRI,"LAGRANGE_MUL_PERIODIC",order); CHKERRQ(ierr);
607  ierr = m_field.set_field_order(0,MBEDGE,"LAGRANGE_MUL_PERIODIC",order); CHKERRQ(ierr);
608  ierr = m_field.set_field_order(0,MBVERTEX,"LAGRANGE_MUL_PERIODIC",1); CHKERRQ(ierr);
609 
610  //build field
611  ierr = m_field.build_fields(); CHKERRQ(ierr);
612 
613  //10 node tets
614  Projection10NodeCoordsOnField ent_method_material(m_field,"MESH_NODE_POSITIONS");
615  ierr = m_field.loop_dofs("MESH_NODE_POSITIONS",ent_method_material,0); CHKERRQ(ierr);
616 
617  //FE
618  ierr = m_field.add_finite_element("PLASTIC"); CHKERRQ(ierr);
619  //Define rows/cols and element data
620  ierr = m_field.modify_finite_element_add_field_row("PLASTIC","DISPLACEMENT"); CHKERRQ(ierr);
621  ierr = m_field.modify_finite_element_add_field_col("PLASTIC","DISPLACEMENT"); CHKERRQ(ierr);
622  ierr = m_field.modify_finite_element_add_field_data("PLASTIC","DISPLACEMENT"); CHKERRQ(ierr);
623  ierr = m_field.modify_finite_element_add_field_data("PLASTIC","MESH_NODE_POSITIONS"); CHKERRQ(ierr);
624 
625  Range newtets;
626  for(_IT_CUBITMESHSETS_BY_BCDATA_TYPE_FOR_LOOP_(m_field,BLOCKSET|UNKNOWNNAME,it)) {
627  if(it->getName() != "MAT_PLASTIC") continue;
628  EntityHandle meshset = it->getMeshset();
629  int id = it->getMeshsetId();
630 
631  rval = m_field.get_moab().get_entities_by_type(
632  meshset,MBTET,newtets,true
633  ); CHKERRQ_MOAB(rval);
634 
635  cout<< "========================== newtets "<<newtets.size()<<endl;
636 
637  //intersection of new and old tets for plastic
638  newtets = intersect(newtets,LatestRefinedTets);
639  ierr = m_field.seed_finite_elements(newtets); CHKERRQ(ierr);
640  }
641 
642  ierr = m_field.add_ents_to_finite_element_by_type(newtets,MBTET,"PLASTIC"); CHKERRQ(ierr);
643  cout<< "========================== newtets "<<newtets.size()<<endl;
644 
645 
646  //================================================================================================================================
647  // Trans-Isotropic Yarns
648  //================================================================================================================================
649 
650  NonlinearElasticElement trans_elastic(m_field,2);
651  trans_elastic.commonData.spatialPositions = "DISPLACEMENT";
652  trans_elastic.commonData.meshPositions = "MESH_NODE_POSITIONS";
653  std::map<int,boost::shared_ptr<SmallStrainTranverslyIsotropicADouble> > tranversly_isotropic_adouble_ptr_map;
654  std::map<int,boost::shared_ptr<SmallStrainTranverslyIsotropicDouble> > tranversly_isotropic_double_ptr_map;
655  bool trans_iso_blocks = false;
656  for(_IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_(m_field,BLOCKSET,it)) {
657  //Get block name
658  string name = it->getName();
659  if (name.compare(0,20,"MAT_ELASTIC_TRANSISO") == 0) {
660  cout<<"================================ it is MAT_ELASTIC_TRANSISO "<<endl;
661  trans_iso_blocks = true;
662  int id = it->getMeshsetId();
663  Mat_Elastic_TransIso mydata;
664  ierr = it->getAttributeDataStructure(mydata); CHKERRQ(ierr);
665  tranversly_isotropic_adouble_ptr_map[id] = boost::make_shared<SmallStrainTranverslyIsotropicADouble>();
666  tranversly_isotropic_double_ptr_map[id] = boost::make_shared<SmallStrainTranverslyIsotropicDouble>();
667  //nu_p, nu_pz, E_p, E_z, G_zp
668  tranversly_isotropic_adouble_ptr_map.at(id)->E_p = mydata.data.Youngp;
669  tranversly_isotropic_double_ptr_map.at(id)->E_p = mydata.data.Youngp;
670  tranversly_isotropic_adouble_ptr_map.at(id)->E_z = mydata.data.Youngz;
671  tranversly_isotropic_double_ptr_map.at(id)->E_z = mydata.data.Youngz;
672  tranversly_isotropic_adouble_ptr_map.at(id)->nu_p = mydata.data.Poissonp;
673  tranversly_isotropic_double_ptr_map.at(id)->nu_p = mydata.data.Poissonp;
674  tranversly_isotropic_adouble_ptr_map.at(id)->nu_pz = mydata.data.Poissonpz;
675  tranversly_isotropic_double_ptr_map.at(id)->nu_pz = mydata.data.Poissonpz;
676  double shear_zp;
677  if(mydata.data.Shearzp!=0) {
678  shear_zp = mydata.data.Shearzp;
679  } else {
680  shear_zp = mydata.data.Youngz/(2*(1+mydata.data.Poissonpz));
681  }
682  tranversly_isotropic_adouble_ptr_map.at(it->getMeshsetId())->G_zp = shear_zp;
683  tranversly_isotropic_double_ptr_map.at(it->getMeshsetId())->G_zp = shear_zp;
684  //get tets from block where material is defined
685  EntityHandle meshset = it->getMeshset();
686  rval = m_field.get_moab().get_entities_by_type(
687  meshset,MBTET,trans_elastic.setOfBlocks[id].tEts,true
688  ); CHKERRQ_MOAB(rval);
689 
690 // cout<<"================================ trans_elastic.setOfBlocks[id].tEts "<<trans_elastic.setOfBlocks[id].tEts.size()<<endl;
691  //intersection of new and old tets for tran-iso part
692  trans_elastic.setOfBlocks[id].tEts = intersect(trans_elastic.setOfBlocks[id].tEts,LatestRefinedTets);
693 // cout<<"================================ trans_elastic.setOfBlocks[id].tEts "<<trans_elastic.setOfBlocks[id].tEts.size()<<endl;
694 
695  //adding material to nonlinear class
696  trans_elastic.setOfBlocks[id].iD = id;
697  //note that material parameters are defined internally in material model
698  trans_elastic.setOfBlocks[id].E = 0; // this is not working for this material
699  trans_elastic.setOfBlocks[id].PoissonRatio = 0; // this is not working for this material
700  trans_elastic.setOfBlocks[id].materialDoublePtr = tranversly_isotropic_double_ptr_map.at(id);
701  trans_elastic.setOfBlocks[id].materialAdoublePtr = tranversly_isotropic_adouble_ptr_map.at(id);
702  ierr = m_field.seed_finite_elements(trans_elastic.setOfBlocks[id].tEts); CHKERRQ(ierr);
703  }
704  }
705 
706 
707  ierr = m_field.add_finite_element("TRAN_ISOTROPIC_ELASTIC"); CHKERRQ(ierr);
708  ierr = m_field.add_finite_element("TRAN_ISOTROPIC_ELASTIC",MF_ZERO); CHKERRQ(ierr);
709  ierr = m_field.modify_finite_element_add_field_row("TRAN_ISOTROPIC_ELASTIC","DISPLACEMENT"); CHKERRQ(ierr);
710  ierr = m_field.modify_finite_element_add_field_col("TRAN_ISOTROPIC_ELASTIC","DISPLACEMENT"); CHKERRQ(ierr);
711  ierr = m_field.modify_finite_element_add_field_data("TRAN_ISOTROPIC_ELASTIC","DISPLACEMENT"); CHKERRQ(ierr);
712  ierr = m_field.modify_finite_element_add_field_data("TRAN_ISOTROPIC_ELASTIC","POTENTIAL_FIELD"); CHKERRQ(ierr);
713 
714  if(m_field.check_field("MESH_NODE_POSITIONS")) {
715  ierr = m_field.modify_finite_element_add_field_data("TRAN_ISOTROPIC_ELASTIC","MESH_NODE_POSITIONS"); CHKERRQ(ierr);
716  }
717  for(
718  map<int,NonlinearElasticElement::BlockData>::iterator sit = trans_elastic.setOfBlocks.begin();
719  sit!=trans_elastic.setOfBlocks.end();sit++
720  ) {
721 
722 
723  cout<<" sit->second.tEts "<<sit->second.tEts.size()<<endl;
724  ierr = m_field.add_ents_to_finite_element_by_type(sit->second.tEts,MBTET,"TRAN_ISOTROPIC_ELASTIC"); CHKERRQ(ierr);
725  }
726 
727 
728  //Rhs
729  trans_elastic.feRhs.getOpPtrVector().push_back(new NonlinearElasticElement::OpGetCommonDataAtGaussPts("DISPLACEMENT",trans_elastic.commonData));
730  trans_elastic.feRhs.getOpPtrVector().push_back(new NonlinearElasticElement::OpGetCommonDataAtGaussPts("POTENTIAL_FIELD",trans_elastic.commonData));
731  if(m_field.check_field("MESH_NODE_POSITIONS")) {
732  trans_elastic.feRhs.getOpPtrVector().push_back(new NonlinearElasticElement::OpGetCommonDataAtGaussPts("MESH_NODE_POSITIONS",trans_elastic.commonData));
733  }
734  map<int,NonlinearElasticElement::BlockData>::iterator sit = trans_elastic.setOfBlocks.begin();
735  for(;sit!=trans_elastic.setOfBlocks.end();sit++) {
736  trans_elastic.feRhs.getOpPtrVector().push_back(new NonlinearElasticElement::OpJacobianPiolaKirchhoffStress("DISPLACEMENT",sit->second,trans_elastic.commonData,2,false,false,true));
737  trans_elastic.feRhs.getOpPtrVector().push_back(new NonlinearElasticElement::OpRhsPiolaKirchhoff("DISPLACEMENT",sit->second,trans_elastic.commonData));
738  }
739 
740  //Lhs
741  trans_elastic.feLhs.getOpPtrVector().push_back(new NonlinearElasticElement::OpGetCommonDataAtGaussPts("DISPLACEMENT",trans_elastic.commonData));
742  trans_elastic.feLhs.getOpPtrVector().push_back(new NonlinearElasticElement::OpGetCommonDataAtGaussPts("POTENTIAL_FIELD",trans_elastic.commonData));
743  if(m_field.check_field("MESH_NODE_POSITIONS")) {
744  trans_elastic.feLhs.getOpPtrVector().push_back(new NonlinearElasticElement::OpGetCommonDataAtGaussPts("MESH_NODE_POSITIONS",trans_elastic.commonData));
745  }
746  sit = trans_elastic.setOfBlocks.begin();
747  for(;sit!=trans_elastic.setOfBlocks.end();sit++) {
748  trans_elastic.feLhs.getOpPtrVector().push_back(new NonlinearElasticElement::OpJacobianPiolaKirchhoffStress("DISPLACEMENT",sit->second,trans_elastic.commonData,2,true,false,true));
749  trans_elastic.feLhs.getOpPtrVector().push_back(new NonlinearElasticElement::OpLhsPiolaKirchhoff_dx("DISPLACEMENT","DISPLACEMENT",sit->second,trans_elastic.commonData));
750  }
751 
752 
753 
754 
755  //================================================================================================================================
756  // Interface Cohesive elemetns
757  //================================================================================================================================
758  //FE Interface
759 
760  ierr = m_field.add_finite_element("INTERFACE"); CHKERRQ(ierr);
761  ierr = m_field.modify_finite_element_add_field_row("INTERFACE","DISPLACEMENT"); CHKERRQ(ierr);
762  ierr = m_field.modify_finite_element_add_field_col("INTERFACE","DISPLACEMENT"); CHKERRQ(ierr);
763  ierr = m_field.modify_finite_element_add_field_data("INTERFACE","DISPLACEMENT"); CHKERRQ(ierr);
764  ierr = m_field.modify_finite_element_add_field_data("INTERFACE","MESH_NODE_POSITIONS"); CHKERRQ(ierr);
765 
766  boost::ptr_vector<CohesiveInterfaceElement::PhysicalEquation> interface_materials;
767 
768  //FIXME this in fact allow only for one type of interface,
769  //problem is Young Moduls in interface mayoung_modulusterial
770 
771 
772  double interface_beta, interface_ft, interface_Gf, interface_h;
773  ierr = PetscOptionsGetReal(0,"-interface_beta",&interface_beta,0); CHKERRQ(ierr);
774  ierr = PetscOptionsGetReal(0,"-interface_ft",&interface_ft,0); CHKERRQ(ierr);
775  ierr = PetscOptionsGetReal(0,"-interface_Gf",&interface_Gf,0); CHKERRQ(ierr);
776  ierr = PetscOptionsGetReal(0,"-interface_h",&interface_h,0); CHKERRQ(ierr);
777 
778 
779  for(_IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_(m_field,BLOCKSET,it)) {
780  cout << endl << *it << endl;
781  //Get block name
782  string name = it->getName();
783  if (name.compare(0,10,"MAT_INTERF") == 0) {
784  Mat_Interf mydata;
785  ierr = it->getAttributeDataStructure(mydata); CHKERRQ(ierr);
786  cout << mydata;
787  interface_materials.push_back(new CohesiveInterfaceElement::PhysicalEquation(m_field));
788 // interface_materials.back().h = 0.134; // we used this because E0=35Gpa and Em=4.7Gpa
789  interface_materials.back().youngModulus = mydata.data.alpha;
790 // interface_materials.back().beta = mydata.data.beta;
791 // interface_materials.back().ft = mydata.data.ft;
792 // interface_materials.back().Gf = mydata.data.Gf;
793 
794  interface_materials.back().h = interface_h;
795  interface_materials.back().beta = interface_beta;
796  interface_materials.back().ft = interface_ft;
797  interface_materials.back().Gf = interface_Gf;
798 
799  EntityHandle meshset = it->getMeshset();
800  Range tris; //All the triangles after splitting
801  rval = moab.get_entities_by_type(meshset,MBTRI,tris,true); CHKERRQ_MOAB(rval);
802 // cout<<"tris.size() "<< tris.size() <<endl;
803  Range ents3d;
804  rval = moab.get_adjacencies(tris,3,false,ents3d,moab::Interface::UNION); CHKERRQ_MOAB(rval);
805  interface_materials.back().pRisms = ents3d.subset_by_type(MBPRISM);
806 // cout<<"unite( interface_prisms, ents3d.subset_by_type(MBPRISM) ) "<<unite( interface_prisms, ents3d.subset_by_type(MBPRISM) )<<endl;
807  }
808  }
809 
810  ierr = m_field.seed_finite_elements(interface_materials.back().pRisms); CHKERRQ(ierr);
811 // cout<<"interface_materials.back().pRisms.size() "<< interface_materials.back().pRisms <<endl;
812  ierr = m_field.add_ents_to_finite_element_by_type(interface_materials.back().pRisms,MBPRISM,"INTERFACE"); CHKERRQ(ierr);
813 
814 
815  // cout<<"young_modulus "<<young_modulus<<endl;
816  { //FIXME
817  boost::ptr_vector<CohesiveInterfaceElement::PhysicalEquation>::iterator pit = interface_materials.begin();
818  for(; pit != interface_materials.end();pit++) {
819  pit->youngModulus = young_modulus; //Young's modulus of bulk elastic material (then for interface E0=E/h)
820  }
821  }
822 
823  CohesiveInterfaceElement cohesive_elements(m_field); //make it golbal, so that its other uses can see it
824  ierr = cohesive_elements.addOps("DISPLACEMENT",interface_materials); CHKERRQ(ierr);
825  //================================================================================================================================
826 
827 
828  BCs_RVELagrange_Periodic lagrangian_element_periodic(m_field);
829 // BCs_RVELagrange_Trac lagrangian_element_trac(m_field);
830 
831  lagrangian_element_periodic.addLagrangiangElement("LAGRANGE_ELEM","DISPLACEMENT","LAGRANGE_MUL_PERIODIC","MESH_NODE_POSITIONS",preriodic_prisms);
832 
833 
834  //Adding element LAGRANGE_ELEM_TRANS
835 // lagrangian_element_trac.addLagrangiangElement("LAGRANGE_ELEM_TRANS","DISPLACEMENT","LAGRANGE_MUL_RIGID_TRANS","MESH_NODE_POSITIONS");
836  ierr = m_field.add_finite_element("LAGRANGE_ELEM_TRANS",MF_ZERO); CHKERRQ(ierr);
837  //============================================================================================================
838  //C row as Lagrange_mul_disp and col as DISPLACEMENT
839  ierr = m_field.modify_finite_element_add_field_row("LAGRANGE_ELEM_TRANS","LAGRANGE_MUL_RIGID_TRANS"); CHKERRQ(ierr);
840  ierr = m_field.modify_finite_element_add_field_col("LAGRANGE_ELEM_TRANS","DISPLACEMENT"); CHKERRQ(ierr);
841  //CT col as Lagrange_mul_disp and row as DISPLACEMENT
842  ierr = m_field.modify_finite_element_add_field_col("LAGRANGE_ELEM_TRANS","LAGRANGE_MUL_RIGID_TRANS"); CHKERRQ(ierr);
843  ierr = m_field.modify_finite_element_add_field_row("LAGRANGE_ELEM_TRANS","DISPLACEMENT"); CHKERRQ(ierr);
844  //data
845  ierr = m_field.modify_finite_element_add_field_data("LAGRANGE_ELEM_TRANS","LAGRANGE_MUL_RIGID_TRANS"); CHKERRQ(ierr);
846  ierr = m_field.modify_finite_element_add_field_data("LAGRANGE_ELEM_TRANS","DISPLACEMENT"); CHKERRQ(ierr);
847  //============================================================================================================
848 
849 
850  Range Tri103;
851  for(_IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_(m_field,SIDESET,it)) {
852  if(it->getName().compare(0,12,"AllBoundSurf") == 0 || it->getMeshsetId() == 103) {
853 
854  //total tri including old + new
855  rval = m_field.get_moab().get_entities_by_type(it->meshset,MBTRI,Tri103,true); CHKERRQ_MOAB(rval);
856  ierr = PetscPrintf(PETSC_COMM_WORLD,"Old triangles on SideSet 103 = %d\n",Tri103.size()); CHKERRQ(ierr);
857  //Only triangles in the new bit-level
858  Tri103 = intersect(AllTris, Tri103);
859  ierr = PetscPrintf(PETSC_COMM_WORLD,"New triangles on SideSet 103 = %d\n",Tri103.size()); CHKERRQ(ierr);
860  ierr = m_field.add_ents_to_finite_element_by_type(Tri103,MBTRI,"LAGRANGE_ELEM_TRANS"); CHKERRQ(ierr);
861  }
862  }
863 
864 // build finite elements
865  ierr = m_field.build_finite_elements(); CHKERRQ(ierr);
866  //build adjacencies
867  ierr = m_field.build_adjacencies(problem_bit_level); CHKERRQ(ierr);
868 
869 
870  DMType dm_name = "PLASTIC_PROB";
871  ierr = DMRegister_MoFEM(dm_name); CHKERRQ(ierr);
872 
873  DM dm;
874  ierr = DMCreate(PETSC_COMM_WORLD,&dm);CHKERRQ(ierr);
875  ierr = DMSetType(dm,dm_name);CHKERRQ(ierr);
876 
877  //set dm datastruture which created mofem datastructures
878  ierr = DMMoFEMCreateMoFEM(dm,&m_field,dm_name,problem_bit_level); CHKERRQ(ierr);
879  ierr = DMSetFromOptions(dm); CHKERRQ(ierr);
880  ierr = DMMoFEMSetIsPartitioned(dm,is_partitioned); CHKERRQ(ierr);
881 
882  //add elements to dm
883  ierr = DMMoFEMAddElement(dm,"PLASTIC"); CHKERRQ(ierr);
884  ierr = DMMoFEMAddElement(dm,"TRAN_ISOTROPIC_ELASTIC"); CHKERRQ(ierr);
885  ierr = DMMoFEMAddElement(dm,"INTERFACE"); CHKERRQ(ierr);
886  ierr = DMMoFEMAddElement(dm,"LAGRANGE_ELEM"); CHKERRQ(ierr);
887  ierr = DMMoFEMAddElement(dm,"LAGRANGE_ELEM_TRANS"); CHKERRQ(ierr);
888  ierr = DMSetUp(dm); CHKERRQ(ierr);
889 
890  //create matrices
891  Vec F,D;
892  Mat Aij;
893  ierr = DMCreateGlobalVector_MoFEM(dm,&D); CHKERRQ(ierr);
894  ierr = VecDuplicate(D,&F); CHKERRQ(ierr);
895  ierr = DMCreateMatrix_MoFEM(dm,&Aij); CHKERRQ(ierr);
896  ierr = VecZeroEntries(D); CHKERRQ(ierr);
897  ierr = VecGhostUpdateBegin(D,INSERT_VALUES,SCATTER_FORWARD); CHKERRQ(ierr);
898  ierr = VecGhostUpdateEnd(D,INSERT_VALUES,SCATTER_FORWARD); CHKERRQ(ierr);
899  ierr = DMoFEMMeshToLocalVector(dm,D,INSERT_VALUES,SCATTER_REVERSE); CHKERRQ(ierr);
900  ierr = MatZeroEntries(Aij); CHKERRQ(ierr);
901 
902  vector<Vec> Fvec(6); //jthis will be used to caluclate the homogenised stiffness matix
903  for(int ii = 0;ii<6;ii++) {
904  ierr = VecDuplicate(D,&Fvec[ii]); CHKERRQ(ierr);
905  ierr = VecZeroEntries(Fvec[ii]); CHKERRQ(ierr);
906  ierr = VecGhostUpdateBegin(Fvec[ii],INSERT_VALUES,SCATTER_FORWARD); CHKERRQ(ierr);
907  ierr = VecGhostUpdateEnd(Fvec[ii],INSERT_VALUES,SCATTER_FORWARD); CHKERRQ(ierr);
908  }
909 
910  //assemble Aij and F
911  SmallStrainPlasticity small_strain_plasticity(m_field);
912  {
913  PetscBool bbar = PETSC_TRUE;
914  ierr = PetscOptionsGetBool(0,"-my_bbar",&bbar,0); CHKERRQ(ierr);
915  small_strain_plasticity.commonData.bBar = bbar;
916  }
917  {
918  small_strain_plasticity.feRhs.getOpPtrVector().push_back(
919  new SmallStrainPlasticity::OpGetCommonDataAtGaussPts(
920  "DISPLACEMENT",small_strain_plasticity.commonData
921  )
922  );
923  small_strain_plasticity.feRhs.getOpPtrVector().push_back(
924  new SmallStrainPlasticity::OpCalculateStress(
925  m_field,"DISPLACEMENT",small_strain_plasticity.commonData,cp,false
926  )
927  );
928  small_strain_plasticity.feRhs.getOpPtrVector().push_back(
929  new SmallStrainPlasticity::OpAssembleRhs(
930  "DISPLACEMENT",small_strain_plasticity.commonData
931  )
932  );
933  small_strain_plasticity.feLhs.getOpPtrVector().push_back(
934  new SmallStrainPlasticity::OpGetCommonDataAtGaussPts(
935  "DISPLACEMENT",small_strain_plasticity.commonData
936  )
937  );
938  small_strain_plasticity.feLhs.getOpPtrVector().push_back(
939  new SmallStrainPlasticity::OpCalculateStress(
940  m_field,"DISPLACEMENT",small_strain_plasticity.commonData,cp,false
941  )
942  );
943  small_strain_plasticity.feLhs.getOpPtrVector().push_back(
944  new SmallStrainPlasticity::OpAssembleLhs(
945  "DISPLACEMENT",small_strain_plasticity.commonData
946  )
947  );
948  small_strain_plasticity.feUpdate.getOpPtrVector().push_back(
949  new SmallStrainPlasticity::OpGetCommonDataAtGaussPts(
950  "DISPLACEMENT",small_strain_plasticity.commonData
951  )
952  );
953  small_strain_plasticity.feUpdate.getOpPtrVector().push_back(
954  new SmallStrainPlasticity::OpCalculateStress(
955  m_field,"DISPLACEMENT",small_strain_plasticity.commonData,cp,false
956  )
957  );
958  small_strain_plasticity.feUpdate.getOpPtrVector().push_back(
959  new SmallStrainPlasticity::OpUpdate(
960  "DISPLACEMENT",small_strain_plasticity.commonData
961  )
962  );
963  ierr = small_strain_plasticity.createInternalVariablesTag(); CHKERRQ(ierr);
964  }
965 
966 
967 
968  lagrangian_element_periodic.setRVEBCsOperatorsNonlinear("DISPLACEMENT","LAGRANGE_MUL_PERIODIC","MESH_NODE_POSITIONS",Aij,Fvec,F,given_strain);
969  BCs_RVELagrange_Trac_Rigid_Trans lagrangian_element_rigid_body_trans(m_field);
970  lagrangian_element_rigid_body_trans.setRVEBCsRigidBodyTranOperators("DISPLACEMENT","LAGRANGE_MUL_RIGID_TRANS",Aij, lagrangian_element_periodic.setOfRVEBC);
971 
972 
973  TimeForceScale time_force_scale("-my_load_history",false);
974  lagrangian_element_periodic.methodsOp.push_back(new TimeForceScale("-my_macro_strian_history",false));
975 
976  //Adding elements to DMSnes
977  //Rhs
978  ierr = DMMoFEMSNESSetFunction(dm,"PLASTIC",&small_strain_plasticity.feRhs,PETSC_NULL,PETSC_NULL); CHKERRQ(ierr);
979  ierr = DMMoFEMSNESSetFunction(dm,"TRAN_ISOTROPIC_ELASTIC",&trans_elastic.feRhs,PETSC_NULL,PETSC_NULL); CHKERRQ(ierr);
980  ierr = DMMoFEMSNESSetFunction(dm,"INTERFACE",&cohesive_elements.feRhs,PETSC_NULL,PETSC_NULL); CHKERRQ(ierr);
981  ierr = DMMoFEMSNESSetFunction(dm,"LAGRANGE_ELEM",&lagrangian_element_periodic.feRVEBCRhs,PETSC_NULL,PETSC_NULL); CHKERRQ(ierr);
982  ierr = DMMoFEMSNESSetFunction(dm,"LAGRANGE_ELEM",&lagrangian_element_periodic.feRVEBCRhsResidual,PETSC_NULL,PETSC_NULL); CHKERRQ(ierr);
983 
984  //Lhs
985  ierr = DMMoFEMSNESSetJacobian(dm,"PLASTIC",&small_strain_plasticity.feLhs,NULL,NULL); CHKERRQ(ierr);
986  ierr = DMMoFEMSNESSetJacobian(dm,"TRAN_ISOTROPIC_ELASTIC",&trans_elastic.feLhs,NULL,NULL); CHKERRQ(ierr);
987  ierr = DMMoFEMSNESSetJacobian(dm,"INTERFACE",&cohesive_elements.feLhs,NULL,NULL); CHKERRQ(ierr);
988  ierr = DMMoFEMSNESSetJacobian(dm,"LAGRANGE_ELEM",&lagrangian_element_periodic.feRVEBCLhs,NULL,NULL); CHKERRQ(ierr);
989  ierr = DMMoFEMSNESSetJacobian(dm,"LAGRANGE_ELEM_TRANS",&lagrangian_element_rigid_body_trans.feRVEBCLhs,NULL,NULL); CHKERRQ(ierr);
990 
991  // Create Time Stepping solver
992  SNES snes;
993  SnesCtx *snes_ctx;
994  SNES snes_one_step_only;
995 
996 
997  ierr = SNESCreate(PETSC_COMM_WORLD,&snes); CHKERRQ(ierr);
998  //ierr = SNESSetDM(snes,dm); CHKERRQ(ierr);
999  ierr = DMMoFEMGetSnesCtx(dm,&snes_ctx); CHKERRQ(ierr);
1000  ierr = SNESSetFunction(snes,F,SnesRhs,snes_ctx); CHKERRQ(ierr);
1001  ierr = SNESSetJacobian(snes,Aij,Aij,SnesMat,snes_ctx); CHKERRQ(ierr);
1002  ierr = SNESSetFromOptions(snes); CHKERRQ(ierr);
1003 
1004  ierr = SNESCreate(PETSC_COMM_WORLD,&snes_one_step_only); CHKERRQ(ierr);
1005  ierr = SNESSetFunction(snes_one_step_only,F,SnesRhs,snes_ctx); CHKERRQ(ierr);
1006  ierr = SNESSetJacobian(snes_one_step_only,Aij,Aij,SnesMat,snes_ctx); CHKERRQ(ierr);
1007  ierr = SNESSetFromOptions(snes_one_step_only); CHKERRQ(ierr);
1008  double atol,rtol,stol;
1009  int maxit,maxf;
1010  SNESGetTolerances(snes_one_step_only,&atol,&rtol,&stol,&maxit,&maxf);
1011  maxit = 1;
1012  atol *= 0;
1013  SNESSetTolerances(snes_one_step_only,atol,rtol,stol,maxit,maxf);
1014  SNESLineSearch linesearch;
1015  SNESGetLineSearch(snes_one_step_only,&linesearch);
1016  SNESLineSearchSetType(linesearch,SNESLINESEARCHBASIC);
1017 
1018 
1019 
1020  PostProcVolumeOnRefinedMesh post_proc(m_field);
1021  ierr = post_proc.generateReferenceElementMesh(); CHKERRQ(ierr);
1022  ierr = post_proc.addFieldValuesPostProc("DISPLACEMENT"); CHKERRQ(ierr);
1023  ierr = post_proc.addFieldValuesGradientPostProc("DISPLACEMENT"); CHKERRQ(ierr);
1024  ierr = post_proc.addFieldValuesPostProc("MESH_NODE_POSITIONS"); CHKERRQ(ierr);
1025 
1026  // Volume calculation
1027  //==============================================================================================================================
1028  VolumeElementForcesAndSourcesCore vol_elem(m_field);
1029  Vec volume_vec;
1030  int volume_vec_ghost[] = { 0 };
1031  ierr = VecCreateGhost(
1032  PETSC_COMM_WORLD,(!m_field.get_comm_rank())?1:0,1,1,volume_vec_ghost,&volume_vec
1033  ); CHKERRQ(ierr);
1034  ierr = VecZeroEntries(volume_vec); CHKERRQ(ierr);
1035  vol_elem.getOpPtrVector().push_back(new VolumeCalculation("DISPLACEMENT",volume_vec));
1036 
1037  ierr = m_field.loop_finite_elements("PLASTIC_PROB","PLASTIC", vol_elem); CHKERRQ(ierr);
1038  ierr = m_field.loop_finite_elements("PLASTIC_PROB","TRAN_ISOTROPIC_ELASTIC", vol_elem); CHKERRQ(ierr);
1039 
1040  ierr = VecAssemblyBegin(volume_vec); CHKERRQ(ierr);
1041  ierr = VecAssemblyEnd(volume_vec); CHKERRQ(ierr);
1042  double rve_volume;
1043  ierr = VecSum(volume_vec,&rve_volume); CHKERRQ(ierr);
1044 
1045  // ierr = VecView(volume_vec,PETSC_VIEWER_STDOUT_WORLD); CHKERRQ(ierr);
1046  ierr = PetscPrintf(PETSC_COMM_WORLD,"RVE Volume %3.2g\n",rve_volume); CHKERRQ(ierr);
1047  ierr = VecDestroy(&volume_vec);
1048  //=============================================================================================================================
1049 
1050  double final_time = 1,delta_time = 0.1;
1051  ierr = PetscOptionsGetReal(0,"-my_final_time",&final_time,0); CHKERRQ(ierr);
1052  ierr = PetscOptionsGetReal(0,"-my_delta_time",&delta_time,0); CHKERRQ(ierr);
1053  double delta_time0 = delta_time;
1054 
1055 // ierr = MatView(Aij,PETSC_VIEWER_DRAW_SELF); CHKERRQ(ierr);
1056 // std::string wait;
1057 // std::cin >> wait;
1058 
1059  Vec D0;
1060  ierr = VecDuplicate(D,&D0); CHKERRQ(ierr);
1061 
1062  int step = 0;
1063  double t = 0;
1064  SNESConvergedReason reason = SNES_CONVERGED_ITERATING;
1065  for(;t<final_time;step++) {
1066  t += delta_time;
1067 
1068 // if(t>final_time){
1069 // break;
1070 // }
1071  PetscPrintf(
1072  PETSC_COMM_WORLD,"Step %d Time %6.4g final time %3.2g\n",step,t,final_time
1073  );
1074  //set time
1075  lagrangian_element_periodic.getLoopFeRVEBCRhs().ts_t = t;
1076  //solve problem at step
1077  ierr = VecAssemblyBegin(D);
1078  ierr = VecAssemblyEnd(D);
1079  ierr = VecCopy(D,D0); CHKERRQ(ierr);
1080  if(step == 0 || reason < 0) {
1081  ierr = SNESSetLagJacobian(snes,-2); CHKERRQ(ierr);
1082  } else {
1083  ierr = SNESSetLagJacobian(snes,-1); CHKERRQ(ierr);
1084  }
1085  ierr = SNESSetLagJacobian(snes_one_step_only,-2);
1086  ierr = SNESSolve(snes_one_step_only,PETSC_NULL,D); CHKERRQ(ierr);
1087  ierr = SNESGetConvergedReason(snes_one_step_only,&reason); CHKERRQ(ierr);
1088  if(reason<0) {
1089  ierr = SNESSolve(snes,PETSC_NULL,D); CHKERRQ(ierr);
1090  }
1091  int its;
1092  ierr = SNESGetIterationNumber(snes,&its); CHKERRQ(ierr);
1093 
1094  ierr = PetscPrintf(
1095  PETSC_COMM_WORLD,"number of Newton iterations = %D\n",its
1096  ); CHKERRQ(ierr);
1097 
1098  ierr = SNESGetConvergedReason(snes,&reason); CHKERRQ(ierr);
1099 
1100  if(reason<0) {
1101  t -= delta_time;
1102  delta_time *= 0.1;
1103  ierr = VecCopy(D0,D); CHKERRQ(ierr);
1104  } else {const int its_d = 60;
1105  const double gamma = 0.5;
1106  const double max_reudction = 1;
1107  const double min_reduction = 1e-1;
1108  double reduction;
1109  reduction = pow((double)its_d/(double)(its+1),gamma);
1110  if(delta_time >= max_reudction*delta_time0 && reduction > 1) {
1111  reduction = 1;
1112  } else if(delta_time <= min_reduction*delta_time0 && reduction < 1) {
1113  reduction = 1;
1114  }
1115 
1116  ierr = PetscPrintf(
1117  PETSC_COMM_WORLD,
1118  "reduction delta_time = %6.4e delta_time = %6.4e\n",
1119  reduction,delta_time
1120  ); CHKERRQ(ierr);
1121  delta_time *= reduction;
1122  if(reduction>1 && delta_time < min_reduction*delta_time0) {
1123  delta_time = min_reduction*delta_time0;
1124  }
1125 
1126  ierr = DMoFEMMeshToGlobalVector(
1127  dm,D,INSERT_VALUES,SCATTER_REVERSE
1128  ); CHKERRQ(ierr);
1129 
1130  ierr = DMoFEMLoopFiniteElements(
1131  dm,"PLASTIC",&small_strain_plasticity.feUpdate
1132  ); CHKERRQ(ierr);
1133 
1134 
1135  ierr = DMoFEMLoopFiniteElements(
1136  dm,"INTERFACE",&cohesive_elements.feHistory
1137  ); CHKERRQ(ierr);
1138 
1139  //Save data on mesh
1140  {
1141  ierr = DMoFEMLoopFiniteElements(
1142  dm,"PLASTIC",&post_proc
1143  ); CHKERRQ(ierr);
1144 
1145  string out_file_name;
1146  {
1147  std::ostringstream stm;
1148  stm << "out_matrix_" << step << ".h5m";
1149  out_file_name = stm.str();
1150  }
1151 
1152  ierr = PetscPrintf(
1153  PETSC_COMM_WORLD,"out file %s\n",out_file_name.c_str()
1154  ); CHKERRQ(ierr);
1155 
1156  rval = post_proc.postProcMesh.write_file(
1157  out_file_name.c_str(),"MOAB","PARALLEL=WRITE_PART"
1158  ); CHKERRQ_MOAB(rval);
1159 
1160 
1161 
1162  ierr = DMoFEMLoopFiniteElements(
1163  dm,"TRAN_ISOTROPIC_ELASTIC",&post_proc
1164  ); CHKERRQ(ierr);
1165  {
1166  std::ostringstream stm;
1167  stm << "out_fibres_" << step << ".h5m";
1168  out_file_name = stm.str();
1169  }
1170  ierr = PetscPrintf(
1171  PETSC_COMM_WORLD,"out file %s\n",out_file_name.c_str()
1172  ); CHKERRQ(ierr);
1173 
1174  rval = post_proc.postProcMesh.write_file(
1175  out_file_name.c_str(),"MOAB","PARALLEL=WRITE_PART"
1176  ); CHKERRQ_MOAB(rval);
1177 
1178  }
1179 
1180  //===================================== Homgenised stress (for given strain) ====================================================
1181  VectorDouble StressHomo;
1182  StressHomo.resize(6);
1183  StressHomo.clear();
1184 
1185  //calculate honmogenised stress
1186  //create a vector for 6 components of homogenized stress
1187  Vec stress_homo;
1188  int stress_homo_ghost[] = { 0,1,2,3,4,5,6 };
1189  ierr = VecCreateGhost(
1190  PETSC_COMM_WORLD,(!m_field.get_comm_rank())?6:0,6,6,stress_homo_ghost,&stress_homo
1191  ); CHKERRQ(ierr);
1192 
1193  lagrangian_element_periodic.setRVEBCsHomoStressOperators("DISPLACEMENT","LAGRANGE_MUL_PERIODIC","MESH_NODE_POSITIONS",stress_homo);
1194 
1195 
1196  PetscScalar *avec;
1197  ierr = VecGetArray(stress_homo,&avec); CHKERRQ(ierr);
1198  ierr = VecZeroEntries(stress_homo); CHKERRQ(ierr);
1199  ierr = m_field.loop_finite_elements(
1200  "PLASTIC_PROB","LAGRANGE_ELEM",lagrangian_element_periodic.getLoopFeRVEBCStress()
1201  ); CHKERRQ(ierr);
1202  ierr = PetscOptionsGetReal(
1203  PETSC_NULL,"-my_rve_volume",&rve_volume,PETSC_NULL
1204  ); CHKERRQ(ierr);
1205  ierr = VecAssemblyBegin(stress_homo); CHKERRQ(ierr);
1206  ierr = VecAssemblyEnd(stress_homo); CHKERRQ(ierr);
1207  ierr = VecScale(stress_homo,1.0/rve_volume); CHKERRQ(ierr);
1208  ierr = VecGhostUpdateBegin(stress_homo,INSERT_VALUES,SCATTER_FORWARD); CHKERRQ(ierr);
1209  ierr = VecGhostUpdateEnd(stress_homo,INSERT_VALUES,SCATTER_FORWARD); CHKERRQ(ierr);
1210 
1211  for(int jj=0; jj<6; jj++) {
1212  StressHomo(jj) = avec[jj];
1213  }
1214 
1215  double scale;
1216  ierr = time_force_scale.getForceScale(t,scale); CHKERRQ(ierr);
1217 
1218  PetscPrintf(PETSC_COMM_WORLD,
1219  "Macro_Strain Homo_Stress Path "
1220  );
1221 
1222  //cout<<"Macro-strain ";
1223  for(int ii=0; ii<6; ii++) {
1224  PetscPrintf(
1225  PETSC_COMM_WORLD,
1226  "%8.5e\t",
1227  t*given_strain(ii)
1228  );
1229  }
1230 
1231  //cout<<"Homo stress ";
1232  for(int ii=0; ii<6; ii++) {
1233  PetscPrintf(
1234  PETSC_COMM_WORLD,
1235  "%8.5e\t",
1236  StressHomo(ii)
1237  );
1238  }
1239 
1240  PetscPrintf(PETSC_COMM_WORLD,
1241  "\n"
1242  );
1243 // //==============================================================================================================================
1244  }
1245  }
1246 
1247  ierr = VecDestroy(&D0); CHKERRQ(ierr);
1248  ierr = MatDestroy(&Aij); CHKERRQ(ierr);
1249  ierr = VecDestroy(&F); CHKERRQ(ierr);
1250  ierr = VecDestroy(&D); CHKERRQ(ierr);
1251 
1252  }
1253  CATCH_ERRORS;
1254 
1255 
1256  MoFEM::Core::Finalize();
1257 }
SmallStrainPlasticity::OpGetCommonDataAtGaussPts
Definition: SmallStrainPlasticity.hpp:413
MethodForForceScaling.hpp
SIDESET
@ SIDESET
Definition: definitions.h:147
MoFEM::CoreInterface::loop_finite_elements
virtual MoFEMErrorCode loop_finite_elements(const std::string problem_name, const std::string &fe_name, FEMethod &method, boost::shared_ptr< NumeredEntFiniteElement_multiIndex > fe_ptr=nullptr, MoFEMTypes bh=MF_EXIST, CacheTupleWeakPtr cache_ptr=CacheTupleSharedPtr(), int verb=DEFAULT_VERBOSITY)=0
Make a loop over finite elements.
MYPCOMM_INDEX
#define MYPCOMM_INDEX
default communicator number PCOMM
Definition: definitions.h:215
SmallStrainParaboloidalPlasticity::nt
double nt
Definition: SmallStrainPlasticityMaterialModels.hpp:175
MoFEM::CoreInterface::loop_dofs
virtual MoFEMErrorCode loop_dofs(const Problem *problem_ptr, const std::string &field_name, RowColData rc, DofMethod &method, int lower_rank, int upper_rank, int verb=DEFAULT_VERBOSITY)=0
Make a loop over dofs.
CohesiveInterfaceElement.hpp
Implementation of linear interface element.
MoFEM::CoreTmp< 0 >
Core (interface) class.
Definition: Core.hpp:82
H1
@ H1
continuous field
Definition: definitions.h:85
BCs_RVELagrange_Periodic
Definition: BCs_RVELagrange_Periodic.hpp:19
SmallStrainPlasticity::feUpdate
MyVolumeFE feUpdate
Definition: SmallStrainPlasticity.hpp:362
SmallStrainParaboloidalPlasticity::lambda
double lambda
Definition: SmallStrainPlasticityMaterialModels.hpp:171
BCs_RVELagrange_Disp.hpp
CohesiveElement::CohesiveInterfaceElement::PhysicalEquation
Constitutive (physical) equation for interface.
Definition: CohesiveInterfaceElement.hpp:51
EntityHandle
Face_CenPos_Handle_multiIndex
multi_index_container< Face_CenPos_Handle, indexed_by< ordered_non_unique< tag< xcoord_tag >, member< Face_CenPos_Handle, double,&Face_CenPos_Handle::xcoord > >, ordered_non_unique< tag< ycoord_tag >, member< Face_CenPos_Handle, double,&Face_CenPos_Handle::ycoord > >, ordered_non_unique< tag< zcoord_tag >, member< Face_CenPos_Handle, double,&Face_CenPos_Handle::zcoord > >, ordered_unique< tag< Tri_Hand_tag >, member< Face_CenPos_Handle, const EntityHandle,&Face_CenPos_Handle::Tri_Hand > >, ordered_unique< tag< Composite_xyzcoord >, composite_key< Face_CenPos_Handle, member< Face_CenPos_Handle, double,&Face_CenPos_Handle::xcoord >, member< Face_CenPos_Handle, double,&Face_CenPos_Handle::ycoord >, member< Face_CenPos_Handle, double,&Face_CenPos_Handle::zcoord > > >> > Face_CenPos_Handle_multiIndex
Definition: homogenisation_periodic_atom_test.cpp:96
young_modulus
double young_modulus
Young modulus.
Definition: plastic.cpp:172
SmallStrainPlasticity::OpUpdate
Definition: SmallStrainPlasticity.hpp:418
PostProcTemplateOnRefineMesh::postProcMesh
moab::Interface & postProcMesh
Definition: PostProcOnRefMesh.hpp:122
BCs_RVELagrange_Periodic::setRVEBCsOperatorsNonlinear
PetscErrorCode setRVEBCsOperatorsNonlinear(string field_name, string lagrang_field_name, string mesh_nodals_positions, Mat aij, vector< Vec > &fvec, Vec f, VectorDouble given_strain)
Definition: BCs_RVELagrange_Periodic.hpp:1091
NonlinearElasticElement::OpJacobianPiolaKirchhoffStress
Operator performs automatic differentiation.
Definition: NonLinearElasticElement.hpp:370
CHKERRQ_MOAB
#define CHKERRQ_MOAB(a)
check error code of MoAB function
Definition: definitions.h:454
MoFEM::Mat_Interf::data
_data_ data
Definition: MaterialBlocks.hpp:441
SmallStrainParaboloidalPlasticity
Small strain plasticity with paraboloidal yield criterion (Isotropic Hardening)
Definition: SmallStrainPlasticityMaterialModels.hpp:162
MoFEM::CoreInterface::modify_finite_element_add_field_row
virtual MoFEMErrorCode modify_finite_element_add_field_row(const std::string &fe_name, const std::string name_row)=0
set field row which finite element use
CohesiveElement::CohesiveInterfaceElement::feLhs
MyPrism feLhs
Definition: CohesiveInterfaceElement.hpp:28
NOBASE
@ NOBASE
Definition: definitions.h:59
SmallStrainPlasticity::ClosestPointProjection::plasticStrain
VectorDouble plasticStrain
Definition: SmallStrainPlasticity.hpp:489
NonlinearElasticElement
structure grouping operators and data used for calculation of nonlinear elastic element
Definition: HookeElement.hpp:27
NonlinearElasticElement::feRhs
MyVolumeFE feRhs
calculate right hand side for tetrahedral elements
Definition: NonLinearElasticElement.hpp:65
_IT_CUBITMESHSETS_BY_BCDATA_TYPE_FOR_LOOP_
#define _IT_CUBITMESHSETS_BY_BCDATA_TYPE_FOR_LOOP_(MESHSET_MANAGER, CUBITBCTYPE, IT)
Iterator that loops over a specific Cubit MeshSet in a moFEM field.
Definition: MeshsetsManager.hpp:49
BCs_RVELagrange_Periodic::methodsOp
boost::ptr_vector< MethodForForceScaling > methodsOp
Definition: BCs_RVELagrange_Periodic.hpp:21
SmallStrainPlasticity::ClosestPointProjection::sTrain
VectorDouble sTrain
Definition: SmallStrainPlasticity.hpp:485
MoFEM::CoreInterface::get_comm_rank
virtual int get_comm_rank() const =0
LAMBDA
#define LAMBDA(E, NU)
Definition: fem_tools.h:22
BCs_RVELagrange_Periodic::feRVEBCLhs
MyPrismFE feRVEBCLhs
Definition: BCs_RVELagrange_Periodic.hpp:34
MoFEM.hpp
CohesiveElement::CohesiveInterfaceElement::feRhs
MyPrism feRhs
Definition: CohesiveInterfaceElement.hpp:27
MoFEM::DMoFEMMeshToLocalVector
PetscErrorCode DMoFEMMeshToLocalVector(DM dm, Vec l, InsertMode mode, ScatterMode scatter_mode)
set local (or ghosted) vector values on mesh for partition only
Definition: DMMoFEM.cpp:527
MoFEM::Projection10NodeCoordsOnField
Projection of edge entities with one mid-node on hierarchical basis.
Definition: Projection10NodeCoordsOnField.hpp:24
MOAB_THROW
#define MOAB_THROW(err)
Check error code of MoAB function and throw MoFEM exception.
Definition: definitions.h:541
MoFEM::CoreTmp< 0 >::Finalize
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition: Core.cpp:112
SmallStrainPlasticity::ClosestPointProjection::tAgs
vector< int > tAgs
Definition: SmallStrainPlasticity.hpp:506
out_file_name
char out_file_name[255]
Definition: initial_diffusion.cpp:53
SmallStrainPlasticity::ClosestPointProjection::snesCreate
PetscErrorCode snesCreate()
Definition: SmallStrainPlasticity.cpp:1027
SmallStrainPlasticity::ClosestPointProjection::createMatAVecR
virtual PetscErrorCode createMatAVecR()
Definition: SmallStrainPlasticity.cpp:924
MoFEM::CoreInterface::add_ents_to_field_by_type
virtual MoFEMErrorCode add_ents_to_field_by_type(const Range &ents, const EntityType type, const std::string &name, int verb=DEFAULT_VERBOSITY)=0
Add entities to field meshset.
MoFEM::DMMoFEMAddElement
PetscErrorCode DMMoFEMAddElement(DM dm, std::string fe_name)
add element to dm
Definition: DMMoFEM.cpp:501
order
constexpr int order
Definition: dg_projection.cpp:18
MoFEM::DeprecatedCoreInterface
Deprecated interface functions.
Definition: DeprecatedCoreInterface.hpp:16
BCs_RVELagrange_Periodic::feRVEBCRhs
MyPrismFE feRVEBCRhs
Definition: BCs_RVELagrange_Periodic.hpp:33
UNKNOWNNAME
@ UNKNOWNNAME
Definition: definitions.h:158
MoFEM::Interface
DeprecatedCoreInterface Interface
Definition: Interface.hpp:1975
BCs_RVELagrange_Periodic::setRVEBCsHomoStressOperators
PetscErrorCode setRVEBCsHomoStressOperators(string field_name, string lagrang_field_name, string mesh_nodals_positions, Vec stress_homo)
Definition: BCs_RVELagrange_Periodic.hpp:1312
MoFEM::DMCreateGlobalVector_MoFEM
PetscErrorCode DMCreateGlobalVector_MoFEM(DM dm, Vec *g)
DMShellSetCreateGlobalVector.
Definition: DMMoFEM.cpp:1171
MoFEM::Mat_Interf
Linear interface data structure.
Definition: MaterialBlocks.hpp:430
MoFEM::CoreInterface::add_ents_to_finite_element_by_type
virtual MoFEMErrorCode add_ents_to_finite_element_by_type(const EntityHandle entities, const EntityType type, const std::string &name, const bool recursive=true)=0
add entities to finite element
Composite_xyzcoord
Definition: homogenisation_periodic_atom_test.cpp:72
MoFEM::Exceptions::rval
static MoFEMErrorCodeGeneric< moab::ErrorCode > rval
Definition: Exceptions.hpp:74
CohesiveElement::CohesiveInterfaceElement::feHistory
MyPrism feHistory
Definition: CohesiveInterfaceElement.hpp:29
SmallStrainParaboloidalPlasticity::sIgma_yt
double sIgma_yt
Definition: SmallStrainPlasticityMaterialModels.hpp:174
MoFEM::DMCreateMatrix_MoFEM
PetscErrorCode DMCreateMatrix_MoFEM(DM dm, Mat *M)
Definition: DMMoFEM.cpp:1201
NonlinearElasticElement::CommonData::spatialPositions
string spatialPositions
Definition: NonLinearElasticElement.hpp:109
xcoord_tag
Definition: homogenisation_periodic_atom_test.cpp:68
NonlinearElasticElement::feLhs
MyVolumeFE feLhs
Definition: NonLinearElasticElement.hpp:67
ContactOps::scale
double scale
Definition: EshelbianContact.hpp:22
SmallStrainPlasticity::ClosestPointProjection::tOl
double tOl
Definition: SmallStrainPlasticity.hpp:491
MoFEM::CoreInterface::add_finite_element
virtual MoFEMErrorCode add_finite_element(const std::string &fe_name, enum MoFEMTypes bh=MF_EXCL, int verb=DEFAULT_VERBOSITY)=0
add finite element
SmallStrainParaboloidalPlasticity::nc
double nc
Definition: SmallStrainPlasticityMaterialModels.hpp:175
MoFEM::CoreInterface::get_moab
virtual moab::Interface & get_moab()=0
MoFEM
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:10
MoFEM::CoreInterface::modify_finite_element_add_field_col
virtual MoFEMErrorCode modify_finite_element_add_field_col(const std::string &fe_name, const std::string name_row)=0
set field col which finite element use
BCs_RVELagrange_Periodic::addLagrangiangElement
PetscErrorCode addLagrangiangElement(const string element_name, const string field_name, const string lagrang_field_name, const string mesh_nodals_positions, Range &periodic_prisms)
Definition: BCs_RVELagrange_Periodic.hpp:52
TimeForceScale
Force scale operator for reading two columns.
Definition: TimeForceScale.hpp:18
MoFEM::PetscOptionsGetReal
PetscErrorCode PetscOptionsGetReal(PetscOptions *, const char pre[], const char name[], PetscReal *dval, PetscBool *set)
Definition: DeprecatedPetsc.hpp:152
MoFEM::CoreInterface::build_finite_elements
virtual MoFEMErrorCode build_finite_elements(int verb=DEFAULT_VERBOSITY)=0
Build finite elements.
main
int main(int argc, char *argv[])
Definition: rve_mech_plas_interface_periodic.cpp:122
CohesiveElement::CohesiveInterfaceElement::addOps
MoFEMErrorCode addOps(const std::string field_name, boost::ptr_vector< CohesiveInterfaceElement::PhysicalEquation > &interfaces)
Driver function settting all operators needed for interface element.
Definition: CohesiveInterfaceElement.hpp:568
MoFEM::Mat_Elastic_TransIso::data
_data_ data
Definition: MaterialBlocks.hpp:382
MoFEM::CoreInterface::add_field
virtual MoFEMErrorCode add_field(const std::string &name, const FieldSpace space, const FieldApproximationBase base, const FieldCoefficientsNumber nb_of_coefficients, const TagType tag_type=MB_TAG_SPARSE, const enum MoFEMTypes bh=MF_EXCL, int verb=DEFAULT_VERBOSITY)=0
Add field.
MoFEM::Mat_Elastic_TransIso
Transverse Isotropic material data structure.
Definition: MaterialBlocks.hpp:371
BCs_RVELagrange_Trac_Rigid_Trans.hpp
CohesiveElement::CohesiveInterfaceElement
Cohesive element implementation.
Definition: CohesiveInterfaceElement.hpp:14
SmallStrainPlasticity::feLhs
MyVolumeFE feLhs
Definition: SmallStrainPlasticity.hpp:361
TimeForceScale.hpp
NonlinearElasticElement::OpLhsPiolaKirchhoff_dx
Definition: NonLinearElasticElement.hpp:556
MoFEM::DMoFEMMeshToGlobalVector
PetscErrorCode DMoFEMMeshToGlobalVector(DM dm, Vec g, InsertMode mode, ScatterMode scatter_mode)
set ghosted vector values on all existing mesh entities
Definition: DMMoFEM.cpp:539
MoFEM::DMMoFEMGetSnesCtx
PetscErrorCode DMMoFEMGetSnesCtx(DM dm, MoFEM::SnesCtx **snes_ctx)
get MoFEM::SnesCtx data structure
Definition: DMMoFEM.cpp:1098
MoFEM::DMRegister_MoFEM
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition: DMMoFEM.cpp:47
SmallStrainPlasticityMaterialModels.hpp
Small Strain plasticity material models.
poisson_ratio
double poisson_ratio
Poisson ratio.
Definition: plastic.cpp:173
SmallStrainPlasticity::commonData
CommonData commonData
Definition: SmallStrainPlasticity.hpp:394
NonlinearElasticElement::CommonData::meshPositions
string meshPositions
Definition: NonLinearElasticElement.hpp:110
NonLinearElasticElement.hpp
Operators and data structures for non-linear elastic analysis.
help
static char help[]
Definition: rve_mech_plas_interface_periodic.cpp:71
Projection10NodeCoordsOnField.hpp
MoFEM::DMMoFEMSNESSetJacobian
PetscErrorCode DMMoFEMSNESSetJacobian(DM dm, const char fe_name[], MoFEM::FEMethod *method, MoFEM::BasicMethod *pre_only, MoFEM::BasicMethod *post_only)
set SNES Jacobian evaluation function
Definition: DMMoFEM.cpp:763
PostProcTemplateVolumeOnRefinedMesh::generateReferenceElementMesh
MoFEMErrorCode generateReferenceElementMesh()
Generate reference mesh on single element.
Definition: PostProcOnRefMesh.hpp:301
PostProcTemplateOnRefineMesh::addFieldValuesGradientPostProc
MoFEMErrorCode addFieldValuesGradientPostProc(const std::string field_name, Vec v=PETSC_NULL)
Add operator to post-process L2 or H1 field gradient.
Definition: PostProcOnRefMesh.hpp:195
PostProcOnRefMesh.hpp
Post-process fields on refined mesh.
MoFEM::CoreInterface::check_field
virtual bool check_field(const std::string &name) const =0
check if field is in database
VolumeCalculation.hpp
Operator can be used with any volume element to calculate sum of volumes of all volumes in the set.
Face_CenPos_Handle
Definition: homogenisation_periodic_atom_test.cpp:61
BCs_RVELagrange_Periodic::getLoopFeRVEBCRhs
MyPrismFE & getLoopFeRVEBCRhs()
Definition: BCs_RVELagrange_Periodic.hpp:39
t
constexpr double t
plate stiffness
Definition: plate.cpp:59
BCs_RVELagrange_Trac_Rigid_Trans::setRVEBCsRigidBodyTranOperators
PetscErrorCode setRVEBCsRigidBodyTranOperators(string field_name, string lagrang_field_name, Mat aij, map< int, RVEBC_Data > setOfRVEBC)
Definition: BCs_RVELagrange_Trac_Rigid_Trans.hpp:109
SmallStrainPlasticity
Small strain finite element implementation.
Definition: SmallStrainPlasticity.hpp:320
MoFEM::DMMoFEMCreateMoFEM
PetscErrorCode DMMoFEMCreateMoFEM(DM dm, MoFEM::Interface *m_field_ptr, const char problem_name[], const MoFEM::BitRefLevel bit_level, const MoFEM::BitRefLevel bit_mask=MoFEM::BitRefLevel().set())
Must be called by user to set MoFEM data structures.
Definition: DMMoFEM.cpp:118
MoFEM::VolumeElementForcesAndSourcesCore
Volume finite element base.
Definition: VolumeElementForcesAndSourcesCore.hpp:26
SmallStrainPlasticity.hpp
Operators and data structures for small strain plasticity.
Face_CenPos_Handle::Face_CenPos_Handle
Face_CenPos_Handle(double _xcoord, double _ycoord, double _zcoord, const EntityHandle _Tri_Hand)
Definition: rve_mech_plas_interface_periodic.cpp:83
BCs_RVELagrange_Disp::setOfRVEBC
map< int, RVEBC_Data > setOfRVEBC
maps side set id with appropriate FluxData
Definition: BCs_RVELagrange_Disp.hpp:56
BCs_RVELagrange_Trac_Rigid_Trans
Definition: BCs_RVELagrange_Trac_Rigid_Trans.hpp:18
mesh_file_name
char mesh_file_name[255]
Definition: mesh_smoothing.cpp:23
SmallStrainPlasticity::ClosestPointProjection::internalVariables
VectorDouble internalVariables
Definition: SmallStrainPlasticity.hpp:486
NonlinearElasticElement::commonData
CommonData commonData
Definition: NonLinearElasticElement.hpp:121
SmallStrainParaboloidalPlasticity::nup
double nup
Definition: SmallStrainPlasticityMaterialModels.hpp:172
MoFEM::SnesRhs
PetscErrorCode SnesRhs(SNES snes, Vec x, Vec f, void *ctx)
This is MoFEM implementation for the right hand side (residual vector) evaluation in SNES solver.
Definition: SnesCtx.cpp:27
Range
MoFEM::PetscOptionsGetRealArray
PetscErrorCode PetscOptionsGetRealArray(PetscOptions *, const char pre[], const char name[], PetscReal dval[], PetscInt *nmax, PetscBool *set)
Definition: DeprecatedPetsc.hpp:192
MoFEM::CoreTmp< 0 >::Initialize
static MoFEMErrorCode Initialize(int *argc, char ***args, const char file[], const char help[])
Initializes the MoFEM database PETSc, MOAB and MPI.
Definition: Core.cpp:72
MF_ZERO
@ MF_ZERO
Definition: definitions.h:98
TimeForceScale::getForceScale
MoFEMErrorCode getForceScale(const double ts_t, double &scale)
Definition: TimeForceScale.hpp:89
SmallStrainParaboloidalPlasticity::mu
double mu
Definition: SmallStrainPlasticityMaterialModels.hpp:170
NonlinearElasticElement::OpGetCommonDataAtGaussPts
Definition: NonLinearElasticElement.hpp:362
CATCH_ERRORS
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:372
NonlinearElasticElement::OpRhsPiolaKirchhoff
Definition: NonLinearElasticElement.hpp:520
MoFEM::DeprecatedCoreInterface::seed_ref_level
DEPRECATED MoFEMErrorCode seed_ref_level(const Range &ents, const BitRefLevel &bit, const bool only_tets=true, int verb=-1)
seed entities in the range and their adjacencies in a particular BitRefLevel
Definition: DeprecatedCoreInterface.cpp:24
PostProcVolumeOnRefinedMesh
Post processing.
Definition: PostProcOnRefMesh.hpp:955
_IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_
#define _IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_(MESHSET_MANAGER, CUBITBCTYPE, IT)
Iterator that loops over a specific Cubit MeshSet having a particular BC meshset in a moFEM field.
Definition: MeshsetsManager.hpp:71
MoFEM::Core
CoreTmp< 0 > Core
Definition: Core.hpp:1094
std
Definition: enable_if.hpp:5
BCs_RVELagrange_Periodic.hpp
BLOCKSET
@ BLOCKSET
Definition: definitions.h:148
SmallTransverselyIsotropic.hpp
NonlinearElasticElement::setOfBlocks
std::map< int, BlockData > setOfBlocks
maps block set id with appropriate BlockData
Definition: NonLinearElasticElement.hpp:100
SmallStrainParaboloidalPlasticity::Hc
double Hc
Definition: SmallStrainPlasticityMaterialModels.hpp:173
SmallStrainPlasticity::OpAssembleRhs
Definition: SmallStrainPlasticity.hpp:438
MoFEM::Exceptions::ierr
static MoFEMErrorCodeGeneric< PetscErrorCode > ierr
Definition: Exceptions.hpp:76
SmallStrainPlasticity::createInternalVariablesTag
PetscErrorCode createInternalVariablesTag()
Definition: SmallStrainPlasticity.cpp:34
AINSWORTH_LEGENDRE_BASE
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:60
MoFEM::ForcesAndSourcesCore::getOpPtrVector
boost::ptr_deque< UserDataOperator > & getOpPtrVector()
Use to push back operator for row operator.
Definition: ForcesAndSourcesCore.hpp:83
EigenMatrix::Vec
const FTensor::Tensor2< T, Dim, Dim > Vec
Definition: MatrixFunction.hpp:66
MoFEM::PetscOptionsGetString
PetscErrorCode PetscOptionsGetString(PetscOptions *, const char pre[], const char name[], char str[], size_t size, PetscBool *set)
Definition: DeprecatedPetsc.hpp:172
SmallStrainParaboloidalPlasticity::Ht
double Ht
Definition: SmallStrainPlasticityMaterialModels.hpp:173
Face_CenPos_Handle_multiIndex
multi_index_container< Face_CenPos_Handle, indexed_by< ordered_non_unique< tag< xcoord_tag >, member< Face_CenPos_Handle, double,&Face_CenPos_Handle::xcoord > >, ordered_non_unique< tag< ycoord_tag >, member< Face_CenPos_Handle, double,&Face_CenPos_Handle::ycoord > >, ordered_non_unique< tag< zcoord_tag >, member< Face_CenPos_Handle, double,&Face_CenPos_Handle::zcoord > >, ordered_unique< tag< Tri_Hand_tag >, member< Face_CenPos_Handle, const EntityHandle,&Face_CenPos_Handle::Tri_Hand > >, ordered_unique< tag< Composite_xyzcoord >, composite_key< Face_CenPos_Handle, member< Face_CenPos_Handle, double,&Face_CenPos_Handle::xcoord >, member< Face_CenPos_Handle, double,&Face_CenPos_Handle::ycoord >, member< Face_CenPos_Handle, double,&Face_CenPos_Handle::zcoord > > >> > Face_CenPos_Handle_multiIndex
Definition: rve_mech_plas_interface_periodic.cpp:115
Tri_Hand_tag
Definition: homogenisation_periodic_atom_test.cpp:71
MoFEM::CoreInterface::build_fields
virtual MoFEMErrorCode build_fields(int verb=DEFAULT_VERBOSITY)=0
MoFEM::CoreInterface::modify_finite_element_add_field_data
virtual MoFEMErrorCode modify_finite_element_add_field_data(const std::string &fe_name, const std::string name_filed)=0
set finite element field data
SmallStrainParaboloidalPlasticity::sIgma_yc
double sIgma_yc
Definition: SmallStrainPlasticityMaterialModels.hpp:174
MoFEM::Types::VectorDouble
UBlasVector< double > VectorDouble
Definition: Types.hpp:68
ReactionDiffusionEquation::D
const double D
diffusivity
Definition: reaction_diffusion.cpp:20
CohesiveElement
Definition: arc_length_interface.cpp:29
SmallStrainPlasticity::CommonData::bBar
bool bBar
Definition: SmallStrainPlasticity.hpp:389
MoFEM::Types::BitRefLevel
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
Definition: Types.hpp:40
SmallStrainPlasticity::OpAssembleLhs
Definition: SmallStrainPlasticity.hpp:456
MoFEM::CoreInterface::build_adjacencies
virtual MoFEMErrorCode build_adjacencies(const Range &ents, int verb=DEFAULT_VERBOSITY)=0
build adjacencies
SmallStrainPlasticity::feRhs
MyVolumeFE feRhs
calculate right hand side for tetrahedral elements
Definition: SmallStrainPlasticity.hpp:360
MoFEM::SnesCtx
Interface for nonlinear (SNES) solver.
Definition: SnesCtx.hpp:13
MoFEM::CoreInterface::set_field_order
virtual MoFEMErrorCode set_field_order(const EntityHandle meshset, const EntityType type, const std::string &name, const ApproximationOrder order, int verb=DEFAULT_VERBOSITY)=0
Set order approximation of the entities in the field.
VolumeCalculation
Calculate volume.
Definition: VolumeCalculation.hpp:17
MoFEM::SnesMat
PetscErrorCode SnesMat(SNES snes, Vec x, Mat A, Mat B, void *ctx)
This is MoFEM implementation for the left hand side (tangent matrix) evaluation in SNES solver.
Definition: SnesCtx.cpp:139
BCs_RVELagrange_Trac.hpp
MoFEM::DMoFEMLoopFiniteElements
PetscErrorCode DMoFEMLoopFiniteElements(DM dm, const char fe_name[], MoFEM::FEMethod *method, CacheTupleWeakPtr cache_ptr=CacheTupleSharedPtr())
Executes FEMethod for finite elements in DM.
Definition: DMMoFEM.cpp:590
MoFEM::PetscOptionsGetInt
PetscErrorCode PetscOptionsGetInt(PetscOptions *, const char pre[], const char name[], PetscInt *ivalue, PetscBool *set)
Definition: DeprecatedPetsc.hpp:142
ycoord_tag
Definition: homogenisation_periodic_atom_test.cpp:69
BCs_RVELagrange_Periodic::getLoopFeRVEBCStress
MyPrismFE & getLoopFeRVEBCStress()
Definition: BCs_RVELagrange_Periodic.hpp:42
BCs_RVELagrange_Disp::feRVEBCLhs
MyTriFE feRVEBCLhs
Definition: BCs_RVELagrange_Disp.hpp:30
MU
#define MU(E, NU)
Definition: fem_tools.h:23
F
@ F
Definition: free_surface.cpp:394
BCs_RVELagrange_Periodic::feRVEBCRhsResidual
MyPrismFE feRVEBCRhsResidual
Definition: BCs_RVELagrange_Periodic.hpp:36
NOFIELD
@ NOFIELD
scalar or vector of scalars describe (no true field)
Definition: definitions.h:84
MoFEM::DMMoFEMSetIsPartitioned
PetscErrorCode DMMoFEMSetIsPartitioned(DM dm, PetscBool is_partitioned)
Definition: DMMoFEM.cpp:1127
MoFEM::DMMoFEMSNESSetFunction
PetscErrorCode DMMoFEMSNESSetFunction(DM dm, const char fe_name[], MoFEM::FEMethod *method, MoFEM::BasicMethod *pre_only, MoFEM::BasicMethod *post_only)
set SNES residual evaluation function
Definition: DMMoFEM.cpp:722
MoFEM::PetscOptionsGetBool
PetscErrorCode PetscOptionsGetBool(PetscOptions *, const char pre[], const char name[], PetscBool *bval, PetscBool *set)
Definition: DeprecatedPetsc.hpp:182
PostProcTemplateOnRefineMesh::addFieldValuesPostProc
MoFEMErrorCode addFieldValuesPostProc(const std::string field_name, Vec v=PETSC_NULL)
Add operator to post-process L2, H1, Hdiv, Hcurl field value.
Definition: PostProcOnRefMesh.hpp:153
zcoord_tag
Definition: homogenisation_periodic_atom_test.cpp:70
SmallStrainPlasticity::OpCalculateStress
Definition: SmallStrainPlasticity.hpp:603
Generated by Doxygen 1.8.17 and hosted at University of Glasgow