v0.14.0
rve_mech_plas_interface_disp.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 <BCs_RVELagrange_Disp.hpp>
39 #include <VolumeCalculation.hpp>
40 
41 #include <CohesiveInterfaceElement.hpp>
42 using namespace CohesiveElement;
43 
44 #ifndef WITH_ADOL_C
45 #error "MoFEM need to be compiled with ADOL-C"
46 #endif
47 #include <adolc/adolc.h>
48 #include <NonLinearElasticElement.hpp>
49 #include <SmallTransverselyIsotropic.hpp>
50 
51 #include <adolc/adolc.h>
52 #include <SmallStrainPlasticity.hpp>
53 #include <SmallStrainPlasticityMaterialModels.hpp>
54 
55 
56 #include <Projection10NodeCoordsOnField.hpp>
57 #include <PostProcOnRefMesh.hpp>
58 
59 #include <string>
60 
61 using namespace boost::numeric;
62 
63 
64 
65 
66 static char help[] =
67 "...\n"
68 "\n";
69 
70 int main(int argc, char *argv[]) {
71  MoFEM::Core::Initialize(&argc,&argv,(char *)0,help);
72 
73  try {
74 
75  moab::Core mb_instance;
76  moab::Interface& moab = mb_instance;
77  ParallelComm* pcomm = ParallelComm::get_pcomm(&moab,MYPCOMM_INDEX);
78  if(pcomm == NULL) pcomm = new ParallelComm(&moab,PETSC_COMM_WORLD);
79 
80  PetscBool flg = PETSC_TRUE;
81  char mesh_file_name[255];
82  ierr = PetscOptionsGetString(PETSC_NULL,"-my_file",mesh_file_name,255,&flg); CHKERRQ(ierr);
83  if(flg != PETSC_TRUE) {
84  SETERRQ(PETSC_COMM_SELF,1,"*** ERROR -my_file (MESH FILE NEEDED)");
85  }
86 
87  PetscInt order;
88  ierr = PetscOptionsGetInt(PETSC_NULL,"-my_order",&order,&flg); CHKERRQ(ierr);
89  if(flg != PETSC_TRUE) {
90  order = 2;
91  }
92  PetscInt bubble_order;
93  ierr = PetscOptionsGetInt(PETSC_NULL,"-my_bubble_order",&bubble_order,&flg); CHKERRQ(ierr);
94  if(flg != PETSC_TRUE) {
95  bubble_order = order;
96  }
97 
98  // use this if your mesh is partitioned and you run code on parts,
99  // you can solve very big problems
100  PetscBool is_partitioned = PETSC_FALSE;
101  ierr = PetscOptionsGetBool(PETSC_NULL,"-my_is_partitioned",&is_partitioned,&flg); CHKERRQ(ierr);
102 
103  //Applied strain on the RVE (vector of length 6) strain=[xx, yy, zz, xy, xz, zy]^T
104  double mygiven_strain[6];
105  int nmax=6;
106  ierr = PetscOptionsGetRealArray(PETSC_NULL,"-my_given_strain",mygiven_strain,&nmax,&flg); CHKERRQ(ierr);
107  VectorDouble given_strain;
108  given_strain.resize(6);
109  cblas_dcopy(6, &mygiven_strain[0], 1, &given_strain(0), 1);
110  cout<<"given_strain ="<<given_strain<<endl;
111 
112  //Read mesh to MOAB
113  if(is_partitioned == PETSC_TRUE) {
114  const char *option;
115  option = "PARALLEL=BCAST_DELETE;PARALLEL_RESOLVE_SHARED_ENTS;PARTITION=PARALLEL_PARTITION;";
116  rval = moab.load_file(mesh_file_name, 0, option); CHKERRQ_MOAB(rval);
117  } else {
118  const char *option;
119  option = "";
120  rval = moab.load_file(mesh_file_name, 0, option); CHKERRQ_MOAB(rval);
121  }
122 
123 
124 // double young_modulus = 1;
125 // SmallStrainParaboloidalPlasticity cp;
126 // {
127 // cp.tAgs.resize(3);
128 // cp.tAgs[0] = 3;
129 // cp.tAgs[1] = 4;
130 // cp.tAgs[2] = 5;
131 // cp.tOl = 1e-12;
132 //
133 // double poisson_ratio = 0.25;
134 // cp.sIgma_yt = 1;
135 // cp.sIgma_yc = 1;
136 //
137 // cp.Ht = 0.1;
138 // cp.Hc = 0.1;
139 //
140 // cp.nup = 0.3;
141 //
142 //
143 // {
144 // ierr = PetscOptionsGetReal(0,"-my_young_modulus",&young_modulus,0); CHKERRQ(ierr);
145 // ierr = PetscOptionsGetReal(0,"-my_poisson_ratio",&poisson_ratio,0); CHKERRQ(ierr);
146 // cp.mu = MU(young_modulus,poisson_ratio);
147 // cp.lambda = LAMBDA(young_modulus,poisson_ratio);
148 // ierr = PetscOptionsGetReal(0,"-my_sigma_yt",&cp.sIgma_yt,0); CHKERRQ(ierr);
149 // ierr = PetscOptionsGetReal(0,"-my_sigma_yc",&cp.sIgma_yc,0); CHKERRQ(ierr);
150 //
151 // ierr = PetscOptionsGetReal(0,"-my_Ht",&cp.Ht,0); CHKERRQ(ierr);
152 // ierr = PetscOptionsGetReal(0,"-my_Hc",&cp.Hc,0); CHKERRQ(ierr);
153 //
154 // ierr = PetscOptionsGetReal(0,"-my_nup",&cp.nup,0); CHKERRQ(ierr);
155 // }
156 //
157 // PetscPrintf(PETSC_COMM_WORLD,"young_modulus = %4.2e \n",young_modulus);
158 // PetscPrintf(PETSC_COMM_WORLD,"poisson_ratio = %4.2e \n",poisson_ratio);
159 //
160 // PetscPrintf(PETSC_COMM_WORLD,"sIgma_yt = %4.2e \n",cp.sIgma_yt);
161 // PetscPrintf(PETSC_COMM_WORLD,"sIgma_yc = %4.2e \n",cp.sIgma_yt);
162 //
163 // PetscPrintf(PETSC_COMM_WORLD,"nup = %4.2e \n",cp.nup);
164 //
165 // cp.sTrain.resize(6,false);
166 // cp.sTrain.clear();
167 // cp.plasticStrain.resize(6,false);
168 // cp.plasticStrain.clear();
169 // cp.internalVariables.resize(2,false);
170 // cp.internalVariables.clear();
171 // cp.createMatAVecR();
172 // cp.snesCreate();
173 // // ierr = SNESSetFromOptions(cp.sNes); CHKERRQ(ierr);
174 // }
175 
176 
177  double young_modulus = 1;
178  SmallStrainJ2Plasticity cp;
179  {
180  cp.tAgs.resize(3);
181  cp.tAgs[0] = 3;
182  cp.tAgs[1] = 4;
183  cp.tAgs[2] = 5;
184  cp.tOl = 1e-12;
185 
186  double poisson_ratio = 0.25;
187  cp.sIgma_y = 1;
188  cp.H = 0.1;
189  cp.K = 0;
190  cp.phi = 1;
191  {
192  ierr = PetscOptionsGetReal(0,"-my_young_modulus",&young_modulus,0); CHKERRQ(ierr);
193  ierr = PetscOptionsGetReal(0,"-my_poisson_ratio",&poisson_ratio,0); CHKERRQ(ierr);
194  cp.mu = MU(young_modulus,poisson_ratio);
195  cp.lambda = LAMBDA(young_modulus,poisson_ratio);
196  ierr = PetscOptionsGetReal(0,"-my_sigma_y",&cp.sIgma_y,0); CHKERRQ(ierr);
197  ierr = PetscOptionsGetReal(0,"-my_H",&cp.H,0); CHKERRQ(ierr);
198  ierr = PetscOptionsGetReal(0,"-my_K",&cp.K,0); CHKERRQ(ierr);
199  ierr = PetscOptionsGetReal(0,"-my_phi",&cp.phi,0); CHKERRQ(ierr);
200  }
201 
202  PetscPrintf(PETSC_COMM_WORLD,"young_modulus = %4.2e \n",young_modulus);
203  PetscPrintf(PETSC_COMM_WORLD,"poisson_ratio = %4.2e \n",poisson_ratio);
204  PetscPrintf(PETSC_COMM_WORLD,"sIgma_y = %4.2e \n",cp.sIgma_y);
205  PetscPrintf(PETSC_COMM_WORLD,"H = %4.2e \n",cp.H);
206  PetscPrintf(PETSC_COMM_WORLD,"K = %4.2e \n",cp.K);
207  PetscPrintf(PETSC_COMM_WORLD,"phi = %4.2e \n",cp.phi);
208 
209 
210  cp.sTrain.resize(6,false);
211  cp.sTrain.clear();
212  cp.plasticStrain.resize(6,false);
213  cp.plasticStrain.clear();
214  cp.internalVariables.resize(7,false);
215  cp.internalVariables.clear();
216  cp.createMatAVecR();
217  cp.snesCreate();
218  // ierr = SNESSetFromOptions(cp.sNes); CHKERRQ(ierr);
219  }
220 
221 
222 
223 
224 
225 
226  //Create MoFEM (Joseph) database
227  MoFEM::Core core(moab);
228  MoFEM::Interface& m_field = core;
229 
230  vector<BitRefLevel> bit_levels;
231  {
232  Tag th_meshset_info;
233  int def_meshset_info[2] = {0,0};
234  rval = moab.tag_get_handle(
235  "MESHSET_INFO",2,MB_TYPE_INTEGER,th_meshset_info,MB_TAG_CREAT|MB_TAG_SPARSE,&def_meshset_info
236  );
237  int meshset_data[2];
238  EntityHandle root = moab.get_root_set();
239  rval = moab.tag_get_data(th_meshset_info,&root,1,meshset_data); CHKERRQ_MOAB(rval);
240  if(meshset_data[0]==0) {
241  meshset_data[0] = 1;
242  rval = moab.tag_set_data(th_meshset_info,&root,1,meshset_data); CHKERRQ_MOAB(rval);
243 
244  }
245  bit_levels.push_back(BitRefLevel().set(meshset_data[0]-1));
246  }
247 
248 
249  EntityHandle out_meshset;
250  rval = moab.create_meshset(MESHSET_SET,out_meshset); CHKERRQ_MOAB(rval);
251  // ierr = m_field.get_problem_finite_elements_entities("POTENTIAL_PROBLEM","POTENTIAL_ELEM",out_meshset); CHKERRQ(ierr);
252  ierr = m_field.get_entities_by_ref_level(bit_levels.back(),BitRefLevel().set(),out_meshset); CHKERRQ(ierr);
253  Range LatestRefinedTets;
254  rval = moab.get_entities_by_type(out_meshset, MBTET,LatestRefinedTets,true); CHKERRQ_MOAB(rval);
255  Range LatestRefinedPrisms;
256  rval = moab.get_entities_by_type(out_meshset, MBPRISM,LatestRefinedPrisms,true); CHKERRQ_MOAB(rval);
257 
258  cout<< "========================== LatestRefinedPrisms "<<LatestRefinedPrisms.size()<<endl;
259  cout<< "========================== LatestRefinedTets "<<LatestRefinedTets.size()<<endl;
260 
261 
262  BitRefLevel problem_bit_level = bit_levels.back();
263  ierr = m_field.seed_ref_level_3D(0,problem_bit_level); CHKERRQ(ierr);
264 
265  // const clock_t begin_time = clock();
266  ierr = m_field.build_fields(); CHKERRQ(ierr);
267  ierr = m_field.build_finite_elements(); CHKERRQ(ierr);
268 
269 
270 
271  Range prims_on_problem_bit_level;
272  ierr = m_field.get_entities_by_type_and_ref_level(
273  problem_bit_level,BitRefLevel().set(),MBPRISM,prims_on_problem_bit_level
274  ); CHKERRQ(ierr);
275  Range tets_on_problem_bit_level;
276  ierr = m_field.get_entities_by_type_and_ref_level(
277  problem_bit_level,BitRefLevel().set(),MBTET,tets_on_problem_bit_level
278  ); CHKERRQ(ierr);
279 
280  //to create meshset from range
281  EntityHandle meshset_prims_on_problem_bit_level;
282  rval = moab.create_meshset(MESHSET_SET,meshset_prims_on_problem_bit_level); CHKERRQ_MOAB(rval);
283  rval = moab.add_entities(meshset_prims_on_problem_bit_level,prims_on_problem_bit_level); CHKERRQ_MOAB(rval);
284  ierr = m_field.seed_ref_level_MESHSET(meshset_prims_on_problem_bit_level,BitRefLevel().set()); CHKERRQ(ierr);
285 
286 
287  //Fields
288  ierr = m_field.add_field("DISPLACEMENT",H1,AINSWORTH_LEGENDRE_BASE,3); CHKERRQ(ierr);
289  ierr = m_field.add_field("MESH_NODE_POSITIONS",H1,AINSWORTH_LEGENDRE_BASE,3,MB_TAG_SPARSE,MF_ZERO); CHKERRQ(ierr);
290  ierr = m_field.add_field("LAGRANGE_MUL_DISP",H1,AINSWORTH_LEGENDRE_BASE,3); CHKERRQ(ierr);
291 
292  //add entitities (by tets) to the field
293  ierr = m_field.add_ents_to_field_by_type(0,MBTET,"DISPLACEMENT"); CHKERRQ(ierr);
294  ierr = m_field.add_ents_to_field_by_type(0,MBTET,"MESH_NODE_POSITIONS"); CHKERRQ(ierr);
295 
296  for(_IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_(m_field,SIDESET,it)) {
297  if(it->getName().compare(0,12,"AllBoundSurf") == 0 || it->getMeshsetId() == 103) {
298  Range tris;
299  rval = m_field.get_moab().get_entities_by_type(it->meshset,MBTRI,tris,true); CHKERRQ_MOAB(rval);
300  ierr = m_field.add_ents_to_field_by_type(tris,MBTRI,"LAGRANGE_MUL_DISP"); CHKERRQ(ierr);
301  }
302  }
303 
304 
305  //set app. order
306  ierr = m_field.set_field_order(0,MBTET,"DISPLACEMENT",bubble_order); CHKERRQ(ierr);
307  ierr = m_field.set_field_order(0,MBTRI,"DISPLACEMENT",order); CHKERRQ(ierr);
308  ierr = m_field.set_field_order(0,MBEDGE,"DISPLACEMENT",order); CHKERRQ(ierr);
309  ierr = m_field.set_field_order(0,MBVERTEX,"DISPLACEMENT",1); CHKERRQ(ierr);
310 
311  ierr = m_field.set_field_order(0,MBTET,"MESH_NODE_POSITIONS",order>1 ? 2 : 1); CHKERRQ(ierr);
312  ierr = m_field.set_field_order(0,MBTRI,"MESH_NODE_POSITIONS",order>1 ? 2 : 1); CHKERRQ(ierr);
313  ierr = m_field.set_field_order(0,MBEDGE,"MESH_NODE_POSITIONS",order>1 ? 2 : 1); CHKERRQ(ierr);
314  ierr = m_field.set_field_order(0,MBVERTEX,"MESH_NODE_POSITIONS",1); CHKERRQ(ierr);
315 
316  ierr = m_field.set_field_order(0,MBTRI,"LAGRANGE_MUL_DISP",order); CHKERRQ(ierr);
317  ierr = m_field.set_field_order(0,MBEDGE,"LAGRANGE_MUL_DISP",order); CHKERRQ(ierr);
318  ierr = m_field.set_field_order(0,MBVERTEX,"LAGRANGE_MUL_DISP",1); CHKERRQ(ierr);
319 
320  //build field
321  ierr = m_field.build_fields(); CHKERRQ(ierr);
322 
323  //10 node tets
324  Projection10NodeCoordsOnField ent_method_material(m_field,"MESH_NODE_POSITIONS");
325  ierr = m_field.loop_dofs("MESH_NODE_POSITIONS",ent_method_material,0); CHKERRQ(ierr);
326 
327  //FE
328  ierr = m_field.add_finite_element("PLASTIC"); CHKERRQ(ierr);
329  //Define rows/cols and element data
330  ierr = m_field.modify_finite_element_add_field_row("PLASTIC","DISPLACEMENT"); CHKERRQ(ierr);
331  ierr = m_field.modify_finite_element_add_field_col("PLASTIC","DISPLACEMENT"); CHKERRQ(ierr);
332  ierr = m_field.modify_finite_element_add_field_data("PLASTIC","DISPLACEMENT"); CHKERRQ(ierr);
333  ierr = m_field.modify_finite_element_add_field_data("PLASTIC","MESH_NODE_POSITIONS"); CHKERRQ(ierr);
334 
335 
336 
337  Range newtets;
338  for(_IT_CUBITMESHSETS_BY_BCDATA_TYPE_FOR_LOOP_(m_field,BLOCKSET|UNKNOWNNAME,it)) {
339  if(it->getName() != "MAT_PLASTIC") continue;
340  EntityHandle meshset = it->getMeshset();
341  int id = it->getMeshsetId();
342 
343  rval = m_field.get_moab().get_entities_by_type(
344  meshset,MBTET,newtets,true
345  ); CHKERRQ_MOAB(rval);
346 
347  //intersection of new and old tets for plastic
348  newtets = intersect(newtets,LatestRefinedTets);
349  ierr = m_field.seed_finite_elements(newtets); CHKERRQ(ierr);
350  }
351 
352  ierr = m_field.add_ents_to_finite_element_by_type(newtets,MBTET,"PLASTIC"); CHKERRQ(ierr);
353 // cout<< "========================== newtets "<<newtets.size()<<endl;
354 
355 
356  //================================================================================================================================
357  // Trans-Isotropic Yarns
358  //================================================================================================================================
359 
360  NonlinearElasticElement trans_elastic(m_field,2);
361  trans_elastic.commonData.spatialPositions = "DISPLACEMENT";
362  trans_elastic.commonData.meshPositions = "MESH_NODE_POSITIONS";
363  std::map<int,boost::shared_ptr<SmallStrainTranverslyIsotropicADouble> > tranversly_isotropic_adouble_ptr_map;
364  std::map<int,boost::shared_ptr<SmallStrainTranverslyIsotropicDouble> > tranversly_isotropic_double_ptr_map;
365  bool trans_iso_blocks = false;
366  for(_IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_(m_field,BLOCKSET,it)) {
367  //Get block name
368  string name = it->getName();
369  if (name.compare(0,20,"MAT_ELASTIC_TRANSISO") == 0) {
370  cout<<"================================ it is MAT_ELASTIC_TRANSISO "<<endl;
371  trans_iso_blocks = true;
372  int id = it->getMeshsetId();
373  Mat_Elastic_TransIso mydata;
374  ierr = it->getAttributeDataStructure(mydata); CHKERRQ(ierr);
375  tranversly_isotropic_adouble_ptr_map[id] = boost::make_shared<SmallStrainTranverslyIsotropicADouble>();
376  tranversly_isotropic_double_ptr_map[id] = boost::make_shared<SmallStrainTranverslyIsotropicDouble>();
377  //nu_p, nu_pz, E_p, E_z, G_zp
378  tranversly_isotropic_adouble_ptr_map.at(id)->E_p = mydata.data.Youngp;
379  tranversly_isotropic_double_ptr_map.at(id)->E_p = mydata.data.Youngp;
380  tranversly_isotropic_adouble_ptr_map.at(id)->E_z = mydata.data.Youngz;
381  tranversly_isotropic_double_ptr_map.at(id)->E_z = mydata.data.Youngz;
382  tranversly_isotropic_adouble_ptr_map.at(id)->nu_p = mydata.data.Poissonp;
383  tranversly_isotropic_double_ptr_map.at(id)->nu_p = mydata.data.Poissonp;
384  tranversly_isotropic_adouble_ptr_map.at(id)->nu_pz = mydata.data.Poissonpz;
385  tranversly_isotropic_double_ptr_map.at(id)->nu_pz = mydata.data.Poissonpz;
386  double shear_zp;
387  if(mydata.data.Shearzp!=0) {
388  shear_zp = mydata.data.Shearzp;
389  } else {
390  shear_zp = mydata.data.Youngz/(2*(1+mydata.data.Poissonpz));
391  }
392  tranversly_isotropic_adouble_ptr_map.at(it->getMeshsetId())->G_zp = shear_zp;
393  tranversly_isotropic_double_ptr_map.at(it->getMeshsetId())->G_zp = shear_zp;
394  //get tets from block where material is defined
395  EntityHandle meshset = it->getMeshset();
396  rval = m_field.get_moab().get_entities_by_type(
397  meshset,MBTET,trans_elastic.setOfBlocks[id].tEts,true
398  ); CHKERRQ_MOAB(rval);
399 
400  //cout<<"================================ sit->second.tEts "<<sit->second.tEts.size()<<endl;
401  //intersection of new and old tets for tran-iso part
402  trans_elastic.setOfBlocks[id].tEts = intersect(trans_elastic.setOfBlocks[id].tEts,LatestRefinedTets);
403  //cout<<"================================ sit->second.tEts "<<sit->second.tEts.size()<<endl;
404 
405  //adding material to nonlinear class
406  trans_elastic.setOfBlocks[id].iD = id;
407  //note that material parameters are defined internally in material model
408  trans_elastic.setOfBlocks[id].E = 0; // this is not working for this material
409  trans_elastic.setOfBlocks[id].PoissonRatio = 0; // this is not working for this material
410  trans_elastic.setOfBlocks[id].materialDoublePtr = tranversly_isotropic_double_ptr_map.at(id);
411  trans_elastic.setOfBlocks[id].materialAdoublePtr = tranversly_isotropic_adouble_ptr_map.at(id);
412  ierr = m_field.seed_finite_elements(trans_elastic.setOfBlocks[id].tEts); CHKERRQ(ierr);
413  }
414  }
415 
416 
417 
418  ierr = m_field.add_finite_element("TRAN_ISOTROPIC_ELASTIC"); CHKERRQ(ierr);
419  ierr = m_field.add_finite_element("TRAN_ISOTROPIC_ELASTIC",MF_ZERO); CHKERRQ(ierr);
420  ierr = m_field.modify_finite_element_add_field_row("TRAN_ISOTROPIC_ELASTIC","DISPLACEMENT"); CHKERRQ(ierr);
421  ierr = m_field.modify_finite_element_add_field_col("TRAN_ISOTROPIC_ELASTIC","DISPLACEMENT"); CHKERRQ(ierr);
422  ierr = m_field.modify_finite_element_add_field_data("TRAN_ISOTROPIC_ELASTIC","DISPLACEMENT"); CHKERRQ(ierr);
423  ierr = m_field.modify_finite_element_add_field_data("TRAN_ISOTROPIC_ELASTIC","POTENTIAL_FIELD"); CHKERRQ(ierr);
424  if(m_field.check_field("MESH_NODE_POSITIONS")) {
425  ierr = m_field.modify_finite_element_add_field_data("TRAN_ISOTROPIC_ELASTIC","MESH_NODE_POSITIONS"); CHKERRQ(ierr);
426  }
427  for(
428  map<int,NonlinearElasticElement::BlockData>::iterator sit = trans_elastic.setOfBlocks.begin();
429  sit!=trans_elastic.setOfBlocks.end();sit++
430  ) {
431  ierr = m_field.add_ents_to_finite_element_by_type(sit->second.tEts,MBTET,"TRAN_ISOTROPIC_ELASTIC"); CHKERRQ(ierr);
432  }
433 
434 
435  //Rhs
436  trans_elastic.feRhs.getOpPtrVector().push_back(new NonlinearElasticElement::OpGetCommonDataAtGaussPts("DISPLACEMENT",trans_elastic.commonData));
437  trans_elastic.feRhs.getOpPtrVector().push_back(new NonlinearElasticElement::OpGetCommonDataAtGaussPts("POTENTIAL_FIELD",trans_elastic.commonData));
438  if(m_field.check_field("MESH_NODE_POSITIONS")) {
439  trans_elastic.feRhs.getOpPtrVector().push_back(new NonlinearElasticElement::OpGetCommonDataAtGaussPts("MESH_NODE_POSITIONS",trans_elastic.commonData));
440  }
441  map<int,NonlinearElasticElement::BlockData>::iterator sit = trans_elastic.setOfBlocks.begin();
442  for(;sit!=trans_elastic.setOfBlocks.end();sit++) {
443  trans_elastic.feRhs.getOpPtrVector().push_back(new NonlinearElasticElement::OpJacobianPiolaKirchhoffStress("DISPLACEMENT",sit->second,trans_elastic.commonData,2,false,false,true));
444  trans_elastic.feRhs.getOpPtrVector().push_back(new NonlinearElasticElement::OpRhsPiolaKirchhoff("DISPLACEMENT",sit->second,trans_elastic.commonData));
445  }
446 
447  //Lhs
448  trans_elastic.feLhs.getOpPtrVector().push_back(new NonlinearElasticElement::OpGetCommonDataAtGaussPts("DISPLACEMENT",trans_elastic.commonData));
449  trans_elastic.feLhs.getOpPtrVector().push_back(new NonlinearElasticElement::OpGetCommonDataAtGaussPts("POTENTIAL_FIELD",trans_elastic.commonData));
450  if(m_field.check_field("MESH_NODE_POSITIONS")) {
451  trans_elastic.feLhs.getOpPtrVector().push_back(new NonlinearElasticElement::OpGetCommonDataAtGaussPts("MESH_NODE_POSITIONS",trans_elastic.commonData));
452  }
453  sit = trans_elastic.setOfBlocks.begin();
454  for(;sit!=trans_elastic.setOfBlocks.end();sit++) {
455  trans_elastic.feLhs.getOpPtrVector().push_back(new NonlinearElasticElement::OpJacobianPiolaKirchhoffStress("DISPLACEMENT",sit->second,trans_elastic.commonData,2,true,false,true));
456  trans_elastic.feLhs.getOpPtrVector().push_back(new NonlinearElasticElement::OpLhsPiolaKirchhoff_dx("DISPLACEMENT","DISPLACEMENT",sit->second,trans_elastic.commonData));
457  }
458 
459  //================================================================================================================================
460  // Interface Cohesive elemetns
461  //================================================================================================================================
462  //FE Interface
463  ierr = m_field.add_finite_element("INTERFACE"); CHKERRQ(ierr);
464  ierr = m_field.modify_finite_element_add_field_row("INTERFACE","DISPLACEMENT"); CHKERRQ(ierr);
465  ierr = m_field.modify_finite_element_add_field_col("INTERFACE","DISPLACEMENT"); CHKERRQ(ierr);
466  ierr = m_field.modify_finite_element_add_field_data("INTERFACE","DISPLACEMENT"); CHKERRQ(ierr);
467  ierr = m_field.modify_finite_element_add_field_data("INTERFACE","MESH_NODE_POSITIONS"); CHKERRQ(ierr);
468 
469  ierr = m_field.add_ents_to_finite_element_by_bit_ref(problem_bit_level,BitRefLevel().set(),"INTERFACE",MBPRISM); CHKERRQ(ierr);
470 
471  boost::ptr_vector<CohesiveInterfaceElement::PhysicalEquation> interface_materials;
472 
473  //FIXME this in fact allow only for one type of interface,
474  //problem is Young Moduls in interface mayoung_modulusterial
475 
476 
477  double interface_beta, interface_ft, interface_Gf, interface_h;
478  ierr = PetscOptionsGetReal(0,"-interface_beta",&interface_beta,0); CHKERRQ(ierr);
479  ierr = PetscOptionsGetReal(0,"-interface_ft",&interface_ft,0); CHKERRQ(ierr);
480  ierr = PetscOptionsGetReal(0,"-interface_Gf",&interface_Gf,0); CHKERRQ(ierr);
481  ierr = PetscOptionsGetReal(0,"-interface_h",&interface_h,0); CHKERRQ(ierr);
482 
483  for(_IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_(m_field,BLOCKSET,it)) {
484  cout << endl << *it << endl;
485  //Get block name
486  string name = it->getName();
487  if (name.compare(0,10,"MAT_INTERF") == 0) {
488  Mat_Interf mydata;
489  ierr = it->getAttributeDataStructure(mydata); CHKERRQ(ierr);
490  cout << mydata;
491  interface_materials.push_back(new CohesiveInterfaceElement::PhysicalEquation(m_field));
492 // interface_materials.back().h = 0.134; // we used this because E0=35Gpa and Em=4.7Gpa
493  interface_materials.back().youngModulus = mydata.data.alpha;
494 // interface_materials.back().beta = mydata.data.beta;
495 // interface_materials.back().ft = mydata.data.ft;
496 // interface_materials.back().Gf = mydata.data.Gf;
497 
498  interface_materials.back().h = interface_h;
499  interface_materials.back().beta = interface_beta;
500  interface_materials.back().ft = interface_ft;
501  interface_materials.back().Gf = interface_Gf;
502 
503  EntityHandle meshset = it->getMeshset();
504  Range tris;
505  rval = moab.get_entities_by_type(meshset,MBTRI,tris,true); CHKERRQ_MOAB(rval);
506  Range ents3d;
507  rval = moab.get_adjacencies(tris,3,false,ents3d,moab::Interface::UNION); CHKERRQ_MOAB(rval);
508  interface_materials.back().pRisms = ents3d.subset_by_type(MBPRISM);
509  }
510  }
511 
512  // cout<<"young_modulus "<<young_modulus<<endl;
513  { //FIXME
514  boost::ptr_vector<CohesiveInterfaceElement::PhysicalEquation>::iterator pit = interface_materials.begin();
515  for(; pit != interface_materials.end();pit++) {
516  pit->youngModulus = young_modulus; //Young's modulus of bulk elastic material (then for interface E0=E/h)
517  }
518  }
519  CohesiveInterfaceElement cohesive_elements(m_field);
520  ierr = cohesive_elements.addOps("DISPLACEMENT",interface_materials); CHKERRQ(ierr);
521  //================================================================================================================================
522 
523  BCs_RVELagrange_Disp lagrangian_element_disp(m_field);
524  lagrangian_element_disp.addLagrangiangElement("LAGRANGE_ELEM","DISPLACEMENT","LAGRANGE_MUL_DISP","MESH_NODE_POSITIONS");
525 
526 
527 
528 
529  //build finite elements
530  ierr = m_field.build_finite_elements(); CHKERRQ(ierr);
531  //build adjacencies
532  ierr = m_field.build_adjacencies(problem_bit_level); CHKERRQ(ierr);
533 
534 
535  DMType dm_name = "PLASTIC_PROB";
536  ierr = DMRegister_MoFEM(dm_name); CHKERRQ(ierr);
537 
538  DM dm;
539  ierr = DMCreate(PETSC_COMM_WORLD,&dm);CHKERRQ(ierr);
540  ierr = DMSetType(dm,dm_name);CHKERRQ(ierr);
541 
542  //set dm datastruture which created mofem datastructures
543  ierr = DMMoFEMCreateMoFEM(dm,&m_field,dm_name,problem_bit_level); CHKERRQ(ierr);
544  ierr = DMSetFromOptions(dm); CHKERRQ(ierr);
545  ierr = DMMoFEMSetIsPartitioned(dm,is_partitioned); CHKERRQ(ierr);
546  //add elements to dm
547  ierr = DMMoFEMAddElement(dm,"PLASTIC"); CHKERRQ(ierr);
548  ierr = DMMoFEMAddElement(dm,"LAGRANGE_ELEM"); CHKERRQ(ierr);
549  ierr = DMMoFEMAddElement(dm,"TRAN_ISOTROPIC_ELASTIC"); CHKERRQ(ierr);
550  ierr = DMMoFEMAddElement(dm,"INTERFACE"); CHKERRQ(ierr);
551 
552  ierr = DMSetUp(dm); CHKERRQ(ierr);
553 
554  //create matrices
555  Vec F,D;
556  Mat Aij;
557  ierr = DMCreateGlobalVector_MoFEM(dm,&D); CHKERRQ(ierr);
558  ierr = VecDuplicate(D,&F); CHKERRQ(ierr);
559  ierr = DMCreateMatrix_MoFEM(dm,&Aij); CHKERRQ(ierr);
560  ierr = VecZeroEntries(D); CHKERRQ(ierr);
561  ierr = VecGhostUpdateBegin(D,INSERT_VALUES,SCATTER_FORWARD); CHKERRQ(ierr);
562  ierr = VecGhostUpdateEnd(D,INSERT_VALUES,SCATTER_FORWARD); CHKERRQ(ierr);
563  ierr = DMoFEMMeshToLocalVector(dm,D,INSERT_VALUES,SCATTER_REVERSE); CHKERRQ(ierr);
564  ierr = MatZeroEntries(Aij); CHKERRQ(ierr);
565 
566  vector<Vec> Fvec(6); //jthis will be used to caluclate the homogenised stiffness matix
567  for(int ii = 0;ii<6;ii++) {
568  ierr = VecDuplicate(D,&Fvec[ii]); CHKERRQ(ierr);
569  ierr = VecZeroEntries(Fvec[ii]); CHKERRQ(ierr);
570  ierr = VecGhostUpdateBegin(Fvec[ii],INSERT_VALUES,SCATTER_FORWARD); CHKERRQ(ierr);
571  ierr = VecGhostUpdateEnd(Fvec[ii],INSERT_VALUES,SCATTER_FORWARD); CHKERRQ(ierr);
572  }
573 
574  //assemble Aij and F
575  SmallStrainPlasticity small_strain_plasticity(m_field);
576  {
577  PetscBool bbar = PETSC_TRUE;
578  ierr = PetscOptionsGetBool(0,"-my_bbar",&bbar,0); CHKERRQ(ierr);
579  small_strain_plasticity.commonData.bBar = bbar;
580  }
581  {
582  small_strain_plasticity.feRhs.getOpPtrVector().push_back(
583  new SmallStrainPlasticity::OpGetCommonDataAtGaussPts(
584  "DISPLACEMENT",small_strain_plasticity.commonData
585  )
586  );
587  small_strain_plasticity.feRhs.getOpPtrVector().push_back(
588  new SmallStrainPlasticity::OpCalculateStress(
589  m_field,"DISPLACEMENT",small_strain_plasticity.commonData,cp
590  )
591  );
592  small_strain_plasticity.feRhs.getOpPtrVector().push_back(
593  new SmallStrainPlasticity::OpAssembleRhs(
594  "DISPLACEMENT",small_strain_plasticity.commonData
595  )
596  );
597  small_strain_plasticity.feLhs.getOpPtrVector().push_back(
598  new SmallStrainPlasticity::OpGetCommonDataAtGaussPts(
599  "DISPLACEMENT",small_strain_plasticity.commonData
600  )
601  );
602  small_strain_plasticity.feLhs.getOpPtrVector().push_back(
603  new SmallStrainPlasticity::OpCalculateStress(
604  m_field,"DISPLACEMENT",small_strain_plasticity.commonData,cp
605  )
606  );
607  small_strain_plasticity.feLhs.getOpPtrVector().push_back(
608  new SmallStrainPlasticity::OpAssembleLhs(
609  "DISPLACEMENT",small_strain_plasticity.commonData
610  )
611  );
612  small_strain_plasticity.feUpdate.getOpPtrVector().push_back(
613  new SmallStrainPlasticity::OpGetCommonDataAtGaussPts(
614  "DISPLACEMENT",small_strain_plasticity.commonData
615  )
616  );
617  small_strain_plasticity.feUpdate.getOpPtrVector().push_back(
618  new SmallStrainPlasticity::OpCalculateStress(
619  m_field,"DISPLACEMENT",small_strain_plasticity.commonData,cp
620  )
621  );
622  small_strain_plasticity.feUpdate.getOpPtrVector().push_back(
623  new SmallStrainPlasticity::OpUpdate(
624  "DISPLACEMENT",small_strain_plasticity.commonData
625  )
626  );
627  ierr = small_strain_plasticity.createInternalVariablesTag(); CHKERRQ(ierr);
628  }
629 
630  lagrangian_element_disp.setRVEBCsOperatorsNonlinear("DISPLACEMENT","LAGRANGE_MUL_DISP","MESH_NODE_POSITIONS",Aij,Fvec,F,given_strain);
631  TimeForceScale time_force_scale("-my_load_history",false);
632  lagrangian_element_disp.methodsOp.push_back(new TimeForceScale("-my_macro_strian_history",false));
633 
634 
635 
636  //Adding elements to DMSnes
637  //Rhs
638  ierr = DMMoFEMSNESSetFunction(dm,"PLASTIC",&small_strain_plasticity.feRhs,PETSC_NULL,PETSC_NULL); CHKERRQ(ierr);
639  ierr = DMMoFEMSNESSetFunction(dm,"TRAN_ISOTROPIC_ELASTIC",&trans_elastic.feRhs,PETSC_NULL,PETSC_NULL); CHKERRQ(ierr);
640  ierr = DMMoFEMSNESSetFunction(dm,"INTERFACE",&cohesive_elements.feRhs,PETSC_NULL,PETSC_NULL); CHKERRQ(ierr);
641  ierr = DMMoFEMSNESSetFunction(dm,"LAGRANGE_ELEM",&lagrangian_element_disp.feRVEBCRhs,PETSC_NULL,PETSC_NULL); CHKERRQ(ierr);
642  ierr = DMMoFEMSNESSetFunction(dm,"LAGRANGE_ELEM",&lagrangian_element_disp.feRVEBCRhsResidual,PETSC_NULL,PETSC_NULL); CHKERRQ(ierr);
643 
644  //Lhs
645  ierr = DMMoFEMSNESSetJacobian(dm,"PLASTIC",&small_strain_plasticity.feLhs,NULL,NULL); CHKERRQ(ierr);
646  ierr = DMMoFEMSNESSetJacobian(dm,"TRAN_ISOTROPIC_ELASTIC",&trans_elastic.feLhs,NULL,NULL); CHKERRQ(ierr);
647  ierr = DMMoFEMSNESSetJacobian(dm,"INTERFACE",&cohesive_elements.feLhs,NULL,NULL); CHKERRQ(ierr);
648  ierr = DMMoFEMSNESSetJacobian(dm,"LAGRANGE_ELEM",&lagrangian_element_disp.feRVEBCLhs,NULL,NULL); CHKERRQ(ierr);
649 
650  // Create Time Stepping solver
651  SNES snes;
652  SnesCtx *snes_ctx;
653 
654  ierr = SNESCreate(PETSC_COMM_WORLD,&snes); CHKERRQ(ierr);
655  //ierr = SNESSetDM(snes,dm); CHKERRQ(ierr);
656  ierr = DMMoFEMGetSnesCtx(dm,&snes_ctx); CHKERRQ(ierr);
657  ierr = SNESSetFunction(snes,F,SnesRhs,snes_ctx); CHKERRQ(ierr);
658  ierr = SNESSetJacobian(snes,Aij,Aij,SnesMat,snes_ctx); CHKERRQ(ierr);
659  ierr = SNESSetFromOptions(snes); CHKERRQ(ierr);
660 
661 
662  PostProcVolumeOnRefinedMesh post_proc(m_field);
663  ierr = post_proc.generateReferenceElementMesh(); CHKERRQ(ierr);
664  ierr = post_proc.addFieldValuesPostProc("DISPLACEMENT"); CHKERRQ(ierr);
665  ierr = post_proc.addFieldValuesGradientPostProc("DISPLACEMENT"); CHKERRQ(ierr);
666  ierr = post_proc.addFieldValuesPostProc("MESH_NODE_POSITIONS"); CHKERRQ(ierr);
667 
668  // Volume calculation
669  //==============================================================================================================================
670  VolumeElementForcesAndSourcesCore vol_elem(m_field);
671  Vec volume_vec;
672  int volume_vec_ghost[] = { 0 };
673  ierr = VecCreateGhost(
674  PETSC_COMM_WORLD,(!m_field.get_comm_rank())?1:0,1,1,volume_vec_ghost,&volume_vec
675  ); CHKERRQ(ierr);
676  ierr = VecZeroEntries(volume_vec); CHKERRQ(ierr);
677  vol_elem.getOpPtrVector().push_back(new VolumeCalculation("DISPLACEMENT",volume_vec));
678 
679  ierr = m_field.loop_finite_elements("PLASTIC_PROB","PLASTIC", vol_elem); CHKERRQ(ierr);
680  ierr = m_field.loop_finite_elements("PLASTIC_PROB","TRAN_ISOTROPIC_ELASTIC", vol_elem); CHKERRQ(ierr);
681 
682  ierr = VecAssemblyBegin(volume_vec); CHKERRQ(ierr);
683  ierr = VecAssemblyEnd(volume_vec); CHKERRQ(ierr);
684  double rve_volume;
685  ierr = VecSum(volume_vec,&rve_volume); CHKERRQ(ierr);
686 
687  // ierr = VecView(volume_vec,PETSC_VIEWER_STDOUT_WORLD); CHKERRQ(ierr);
688  ierr = PetscPrintf(PETSC_COMM_WORLD,"RVE Volume %3.2g\n",rve_volume); CHKERRQ(ierr);
689  ierr = VecDestroy(&volume_vec);
690  //=============================================================================================================================
691 
692  double final_time = 1,delta_time = 0.1;
693  ierr = PetscOptionsGetReal(0,"-my_final_time",&final_time,0); CHKERRQ(ierr);
694  ierr = PetscOptionsGetReal(0,"-my_delta_time",&delta_time,0); CHKERRQ(ierr);
695  double delta_time0 = delta_time;
696 
697 // ierr = MatView(Aij,PETSC_VIEWER_DRAW_SELF); CHKERRQ(ierr);
698 // std::string wait;
699 // std::cin >> wait;
700 
701  Vec D0;
702  ierr = VecDuplicate(D,&D0); CHKERRQ(ierr);
703 
704  int step = 0;
705  double t = 0;
706  SNESConvergedReason reason = SNES_CONVERGED_ITERATING;
707  for(;t<final_time;step++) {
708  t += delta_time;
709 
710  if(t>final_time){
711  break;
712  }
713  PetscPrintf(
714  PETSC_COMM_WORLD,"Step %d Time %6.4g final time %3.2g\n",step,t,final_time
715  );
716  //set time
717  lagrangian_element_disp.getLoopFeRVEBCRhs().ts_t = t;
718  //solve problem at step
719  ierr = VecAssemblyBegin(D);
720  ierr = VecAssemblyEnd(D);
721  ierr = VecCopy(D,D0); CHKERRQ(ierr);
722  if(step == 0 || reason < 0) {
723  ierr = SNESSetLagJacobian(snes,-2); CHKERRQ(ierr);
724  } else {
725  ierr = SNESSetLagJacobian(snes,-1); CHKERRQ(ierr);
726  }
727  ierr = SNESSolve(snes,PETSC_NULL,D); CHKERRQ(ierr);
728 
729  int its;
730  ierr = SNESGetIterationNumber(snes,&its); CHKERRQ(ierr);
731 
732  ierr = PetscPrintf(
733  PETSC_COMM_WORLD,"number of Newton iterations = %D\n",its
734  ); CHKERRQ(ierr);
735 
736  ierr = SNESGetConvergedReason(snes,&reason); CHKERRQ(ierr);
737 
738  if(reason<0) {
739  t -= delta_time;
740  delta_time *= 0.1;
741  ierr = VecCopy(D0,D); CHKERRQ(ierr);
742  } else {const int its_d = 60;
743  const double gamma = 0.5;
744  const double max_reudction = 1;
745  const double min_reduction = 1e-1;
746  double reduction;
747  reduction = pow((double)its_d/(double)(its+1),gamma);
748  if(delta_time >= max_reudction*delta_time0 && reduction > 1) {
749  reduction = 1;
750  } else if(delta_time <= min_reduction*delta_time0 && reduction < 1) {
751  reduction = 1;
752  }
753 
754  ierr = PetscPrintf(
755  PETSC_COMM_WORLD,
756  "reduction delta_time = %6.4e delta_time = %6.4e\n",
757  reduction,delta_time
758  ); CHKERRQ(ierr);
759  delta_time *= reduction;
760  if(reduction>1 && delta_time < min_reduction*delta_time0) {
761  delta_time = min_reduction*delta_time0;
762  }
763 
764  ierr = DMoFEMMeshToGlobalVector(
765  dm,D,INSERT_VALUES,SCATTER_REVERSE
766  ); CHKERRQ(ierr);
767 
768  ierr = DMoFEMLoopFiniteElements(
769  dm,"PLASTIC",&small_strain_plasticity.feUpdate
770  ); CHKERRQ(ierr);
771 
772 
773 
774  ierr = DMoFEMLoopFiniteElements(
775  dm,"INTERFACE",&cohesive_elements.feHistory
776  ); CHKERRQ(ierr);
777 
778 
779  //Save data on mesh
780  {
781  ierr = DMoFEMLoopFiniteElements(
782  dm,"PLASTIC",&post_proc
783  ); CHKERRQ(ierr);
784  string out_file_name;
785  {
786  std::ostringstream stm;
787  stm << "out_" << step << ".h5m";
788  out_file_name = stm.str();
789  }
790  ierr = PetscPrintf(
791  PETSC_COMM_WORLD,"out file %s\n",out_file_name.c_str()
792  ); CHKERRQ(ierr);
793  rval = post_proc.postProcMesh.write_file(
794  out_file_name.c_str(),"MOAB","PARALLEL=WRITE_PART"
795  ); CHKERRQ_MOAB(rval);
796  }
797 
798  //===================================== Homgenised stress (for given strain) ====================================================
799  VectorDouble StressHomo;
800  StressHomo.resize(6);
801  StressHomo.clear();
802 
803  //calculate honmogenised stress
804  //create a vector for 6 components of homogenized stress
805  Vec stress_homo;
806  int stress_homo_ghost[] = { 0,1,2,3,4,5,6 };
807  ierr = VecCreateGhost(
808  PETSC_COMM_WORLD,(!m_field.get_comm_rank())?6:0,6,6,stress_homo_ghost,&stress_homo
809  ); CHKERRQ(ierr);
810 
811  lagrangian_element_disp.setRVEBCsHomoStressOperators("DISPLACEMENT","LAGRANGE_MUL_DISP","MESH_NODE_POSITIONS",stress_homo);
812 
813 
814  PetscScalar *avec;
815  ierr = VecGetArray(stress_homo,&avec); CHKERRQ(ierr);
816  ierr = VecZeroEntries(stress_homo); CHKERRQ(ierr);
817  ierr = m_field.loop_finite_elements(
818  "PLASTIC_PROB","LAGRANGE_ELEM",lagrangian_element_disp.getLoopFeRVEBCStress()
819  ); CHKERRQ(ierr);
820  ierr = PetscOptionsGetReal(
821  PETSC_NULL,"-my_rve_volume",&rve_volume,PETSC_NULL
822  ); CHKERRQ(ierr);
823  ierr = VecAssemblyBegin(stress_homo); CHKERRQ(ierr);
824  ierr = VecAssemblyEnd(stress_homo); CHKERRQ(ierr);
825  ierr = VecScale(stress_homo,1.0/rve_volume); CHKERRQ(ierr);
826  ierr = VecGhostUpdateBegin(stress_homo,INSERT_VALUES,SCATTER_FORWARD); CHKERRQ(ierr);
827  ierr = VecGhostUpdateEnd(stress_homo,INSERT_VALUES,SCATTER_FORWARD); CHKERRQ(ierr);
828 
829  for(int jj=0; jj<6; jj++) {
830  StressHomo(jj) = avec[jj];
831  }
832 
833  double scale;
834  ierr = time_force_scale.getForceScale(t,scale); CHKERRQ(ierr);
835 
836  PetscPrintf(PETSC_COMM_WORLD,
837  "Macro_Strain Homo_Stress Path "
838  );
839 
840  //cout<<"Macro-strain ";
841  for(int ii=0; ii<6; ii++) {
842  PetscPrintf(
843  PETSC_COMM_WORLD,
844  "%8.5e\t",
845  t*given_strain(ii)
846  );
847  }
848 
849  //cout<<"Homo stress ";
850  for(int ii=0; ii<6; ii++) {
851  PetscPrintf(
852  PETSC_COMM_WORLD,
853  "%8.5e\t",
854  StressHomo(ii)
855  );
856  }
857 
858  PetscPrintf(PETSC_COMM_WORLD,
859  "\n"
860  );
861 // //==============================================================================================================================
862  }
863  }
864  ierr = VecDestroy(&D0); CHKERRQ(ierr);
865 
866  ierr = MatDestroy(&Aij); CHKERRQ(ierr);
867  ierr = VecDestroy(&F); CHKERRQ(ierr);
868  ierr = VecDestroy(&D); CHKERRQ(ierr);
869  }
870  CATCH_ERRORS;
871 
872  MoFEM::Core::Finalize();
873 }
BCs_RVELagrange_Disp::feRVEBCRhs
MyTriFE feRVEBCRhs
Definition: BCs_RVELagrange_Disp.hpp:29
SIDESET
@ SIDESET
Definition: definitions.h:160
BCs_RVELagrange_Disp::getLoopFeRVEBCStress
MyTriFE & getLoopFeRVEBCStress()
Definition: BCs_RVELagrange_Disp.hpp:39
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:228
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_Disp.hpp
CohesiveElement::CohesiveInterfaceElement::PhysicalEquation
Constitutive (physical) equation for interface.
Definition: CohesiveInterfaceElement.hpp:51
EntityHandle
young_modulus
double young_modulus
Young modulus.
Definition: plastic.cpp:121
PostProcTemplateOnRefineMesh::postProcMesh
moab::Interface & postProcMesh
Definition: PostProcOnRefMesh.hpp:122
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:467
MoFEM::Mat_Interf::data
_data_ data
Definition: MaterialBlocks.hpp:441
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
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
help
static char help[]
Definition: rve_mech_plas_interface_disp.cpp:66
MoFEM::CoreInterface::get_comm_rank
virtual int get_comm_rank() const =0
LAMBDA
#define LAMBDA(E, NU)
Definition: fem_tools.h:22
BCs_RVELagrange_Disp::getLoopFeRVEBCRhs
MyTriFE & getLoopFeRVEBCRhs()
Definition: BCs_RVELagrange_Disp.hpp:36
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:523
MoFEM::Projection10NodeCoordsOnField
Projection of edge entities with one mid-node on hierarchical basis.
Definition: Projection10NodeCoordsOnField.hpp:24
MoFEM::CoreTmp< 0 >::Finalize
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition: Core.cpp:112
out_file_name
char out_file_name[255]
Definition: initial_diffusion.cpp:53
scale
double scale
Definition: plastic.cpp:119
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:497
order
constexpr int order
Definition: dg_projection.cpp:18
MoFEM::DeprecatedCoreInterface
Deprecated interface functions.
Definition: DeprecatedCoreInterface.hpp:16
UNKNOWNNAME
@ UNKNOWNNAME
Definition: definitions.h:171
MoFEM::Interface
DeprecatedCoreInterface Interface
Definition: Interface.hpp:2010
MoFEM::DMCreateGlobalVector_MoFEM
PetscErrorCode DMCreateGlobalVector_MoFEM(DM dm, Vec *g)
DMShellSetCreateGlobalVector.
Definition: DMMoFEM.cpp:1167
MoFEM::Mat_Interf
Linear interface data structure.
Definition: MaterialBlocks.hpp:430
BCs_RVELagrange_Disp
Definition: BCs_RVELagrange_Disp.hpp:18
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
MoFEM::Exceptions::rval
static MoFEMErrorCodeGeneric< moab::ErrorCode > rval
Definition: Exceptions.hpp:74
CohesiveElement::CohesiveInterfaceElement::feHistory
MyPrism feHistory
Definition: CohesiveInterfaceElement.hpp:29
BCs_RVELagrange_Disp::setRVEBCsHomoStressOperators
PetscErrorCode setRVEBCsHomoStressOperators(string field_name, string lagrang_field_name, string mesh_nodals_positions, Vec Stress_Homo)
Definition: BCs_RVELagrange_Disp.hpp:671
MoFEM::DMCreateMatrix_MoFEM
PetscErrorCode DMCreateMatrix_MoFEM(DM dm, Mat *M)
Definition: DMMoFEM.cpp:1197
NonlinearElasticElement::CommonData::spatialPositions
string spatialPositions
Definition: NonLinearElasticElement.hpp:109
NonlinearElasticElement::feLhs
MyVolumeFE feLhs
Definition: NonLinearElasticElement.hpp:67
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
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
main
int main(int argc, char *argv[])
Definition: rve_mech_plas_interface_disp.cpp:70
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.
BCs_RVELagrange_Disp::methodsOp
boost::ptr_vector< MethodForForceScaling > methodsOp
Definition: BCs_RVELagrange_Disp.hpp:25
MoFEM::CoreInterface::add_ents_to_finite_element_by_bit_ref
virtual MoFEMErrorCode add_ents_to_finite_element_by_bit_ref(const BitRefLevel &bit, const BitRefLevel &mask, const std::string &name, EntityType type, int verb=DEFAULT_VERBOSITY)=0
add TET entities from given refinement level to finite element database given by name
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::Mat_Elastic_TransIso
Transverse Isotropic material data structure.
Definition: MaterialBlocks.hpp:371
CohesiveElement::CohesiveInterfaceElement
Cohesive element implementation.
Definition: CohesiveInterfaceElement.hpp:14
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:535
MoFEM::DMMoFEMGetSnesCtx
PetscErrorCode DMMoFEMGetSnesCtx(DM dm, MoFEM::SnesCtx **snes_ctx)
get MoFEM::SnesCtx data structure
Definition: DMMoFEM.cpp:1094
MoFEM::DMRegister_MoFEM
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition: DMMoFEM.cpp:43
poisson_ratio
double poisson_ratio
Poisson ratio.
Definition: plastic.cpp:122
NonlinearElasticElement::CommonData::meshPositions
string meshPositions
Definition: NonLinearElasticElement.hpp:110
NonLinearElasticElement.hpp
Operators and data structures for non-linear elastic analysis.
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:759
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.
t
constexpr double t
plate stiffness
Definition: plate.cpp:58
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:114
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_field)=0
set finite element field data
MoFEM::VolumeElementForcesAndSourcesCore
Volume finite element base.
Definition: VolumeElementForcesAndSourcesCore.hpp:26
mesh_file_name
char mesh_file_name[255]
Definition: mesh_smoothing.cpp:22
NonlinearElasticElement::commonData
CommonData commonData
Definition: NonLinearElasticElement.hpp:121
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:111
TimeForceScale::getForceScale
MoFEMErrorCode getForceScale(const double ts_t, double &scale)
Definition: TimeForceScale.hpp:89
NonlinearElasticElement::OpGetCommonDataAtGaussPts
Definition: NonLinearElasticElement.hpp:362
CATCH_ERRORS
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:385
NonlinearElasticElement::OpRhsPiolaKirchhoff
Definition: NonLinearElasticElement.hpp:520
PostProcVolumeOnRefinedMesh
Post processing.
Definition: PostProcOnRefMesh.hpp:955
BCs_RVELagrange_Disp::feRVEBCRhsResidual
MyTriFE feRVEBCRhsResidual
Definition: BCs_RVELagrange_Disp.hpp:32
_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:1148
std
Definition: enable_if.hpp:5
BLOCKSET
@ BLOCKSET
Definition: definitions.h:161
SmallTransverselyIsotropic.hpp
NonlinearElasticElement::setOfBlocks
std::map< int, BlockData > setOfBlocks
maps block set id with appropriate BlockData
Definition: NonLinearElasticElement.hpp:100
MoFEM::Exceptions::ierr
static MoFEMErrorCodeGeneric< PetscErrorCode > ierr
Definition: Exceptions.hpp:76
AINSWORTH_LEGENDRE_BASE
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:60
BCs_RVELagrange_Disp::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_Disp.hpp:535
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
MoFEM::CoreInterface::build_fields
virtual MoFEMErrorCode build_fields(int verb=DEFAULT_VERBOSITY)=0
MoFEM::DeprecatedCoreInterface::seed_ref_level_3D
DEPRECATED MoFEMErrorCode seed_ref_level_3D(const EntityHandle meshset, const BitRefLevel &bit, int verb=-1)
seed 2D entities in the meshset and their adjacencies (only TETs adjacencies) in a particular BitRefL...
Definition: DeprecatedCoreInterface.cpp:18
MoFEM::Types::VectorDouble
UBlasVector< double > VectorDouble
Definition: Types.hpp:68
BCs_RVELagrange_Disp::addLagrangiangElement
PetscErrorCode addLagrangiangElement(const string element_name, const string field_name, const string lagrang_field_name, const string mesh_nodals_positions)
Definition: BCs_RVELagrange_Disp.hpp:58
ReactionDiffusionEquation::D
const double D
diffusivity
Definition: reaction_diffusion.cpp:20
CohesiveElement
Definition: arc_length_interface.cpp:29
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
MoFEM::CoreInterface::build_adjacencies
virtual MoFEMErrorCode build_adjacencies(const Range &ents, int verb=DEFAULT_VERBOSITY)=0
build adjacencies
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
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:586
MoFEM::PetscOptionsGetInt
PetscErrorCode PetscOptionsGetInt(PetscOptions *, const char pre[], const char name[], PetscInt *ivalue, PetscBool *set)
Definition: DeprecatedPetsc.hpp:142
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.
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
MoFEM::DMMoFEMSetIsPartitioned
PetscErrorCode DMMoFEMSetIsPartitioned(DM dm, PetscBool is_partitioned)
Definition: DMMoFEM.cpp:1123
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:718
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
Generated by Doxygen 1.8.17 and hosted at University of Glasgow