v0.14.0
homogenisation_trac_atom_test.cpp
Go to the documentation of this file.
1 /* This file is part of MoFEM.
2  * MoFEM is free software: you can redistribute it and/or modify it under
3  * the terms of the GNU Lesser General Public License as published by the
4  * Free Software Foundation, either version 3 of the License, or (at your
5  * option) any later version.
6  *
7  * MoFEM is distributed in the hope that it will be useful, but WITHOUT
8  * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
9  * FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
10  * License for more details.
11  *
12  * You should have received a copy of the GNU Lesser General Public
13  * License along with MoFEM. If not, see <http://www.gnu.org/licenses/>. */
14 
15 #include <MoFEM.hpp>
16 using namespace MoFEM;
17 
18 #include <Projection10NodeCoordsOnField.hpp>
19 #include <petsctime.h>
20 
21 // #include <FEMethod_LowLevelStudent.hpp>
22 // #include <FEMethod_UpLevelStudent.hpp>
23 
24 #include <adolc/adolc.h>
25 #include <NonLinearElasticElement.hpp>
26 #include <Hooke.hpp>
27 
28 #include <boost/shared_ptr.hpp>
29 #include <boost/numeric/ublas/vector_proxy.hpp>
30 #include <boost/iostreams/tee.hpp>
31 #include <boost/iostreams/stream.hpp>
32 #include <fstream>
33 #include <iostream>
34 
35 #include <MethodForForceScaling.hpp>
36 #include <TimeForceScale.hpp>
37 #include <VolumeCalculation.hpp>
38 #include <BCs_RVELagrange_Disp.hpp>
39 #include <BCs_RVELagrange_Trac.hpp>
40 #include <BCs_RVELagrange_Trac_Rigid_Rot.hpp>
41 #include <BCs_RVELagrange_Trac_Rigid_Trans.hpp>
42 
43 
44 #include <boost/ptr_container/ptr_map.hpp>
45 #include <PostProcOnRefMesh.hpp>
46 #include <PostProcStresses.hpp>
47 
48 
49 
50 
51 static char help[] = "...\n\n";
52 #define RND_EPS 1e-6
53 
54 //Rounding
55 double roundn(double n)
56 {
57  //break n into fractional part (fract) and integral part (intp)
58  double fract, intp;
59  fract = modf(n,&intp);
60 
61 // //round up
62 // if (fract>=.5)
63 // {
64 // n*=10;
65 // ceil(n);
66 // n/=10;
67 // }
68 // //round down
69 // if (fract<=.5)
70 // {
71 // n*=10;
72 // floor(n);
73 // n/=10;
74 // }
75  // case where n approximates zero, set n to "positive" zero
76  if (abs(intp)==0)
77  {
78  if(abs(fract)<=RND_EPS)
79  {
80  n=0.000;
81  }
82  }
83  return n;
84 }
85 
86 int main(int argc, char *argv[]) {
87 
88  MoFEM::Core::Initialize(&argc,&argv,(char *)0,help);
89 
90  try {
91 
92  moab::Core mb_instance;
93  moab::Interface& moab = mb_instance;
94  int rank;
95  MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
96 
97  //Reade parameters from line command
98  PetscBool flg = PETSC_TRUE;
99  char mesh_file_name[255];
100  ierr = PetscOptionsGetString(PETSC_NULL,"-my_file",mesh_file_name,255,&flg); CHKERRQ(ierr);
101  if(flg != PETSC_TRUE) {
102  SETERRQ(PETSC_COMM_SELF,1,"*** ERROR -my_file (MESH FILE NEEDED)");
103  }
104  PetscInt order;
105  ierr = PetscOptionsGetInt(PETSC_NULL,"-my_order",&order,&flg); CHKERRQ(ierr);
106  if(flg != PETSC_TRUE) {
107  order = 5;
108  }
109 
110  //Read mesh to MOAB
111  const char *option;
112  option = "";//"PARALLEL=BCAST;";//;DEBUG_IO";
113  rval = moab.load_file(mesh_file_name, 0, option); CHKERRQ_MOAB(rval);
114  ParallelComm* pcomm = ParallelComm::get_pcomm(&moab,MYPCOMM_INDEX);
115  if(pcomm == NULL) pcomm = new ParallelComm(&moab,PETSC_COMM_WORLD);
116 
117  //We need that for code profiling
118  PetscLogDouble t1,t2;
119  PetscLogDouble v1,v2;
120  ierr = PetscTime(&v1); CHKERRQ(ierr);
121  ierr = PetscGetCPUTime(&t1); CHKERRQ(ierr);
122 
123  //Create MoFEM (Joseph) database
124  MoFEM::Core core(moab);
125  MoFEM::Interface& m_field = core;
126 
127  //ref meshset ref level 0
128  ierr = m_field.seed_ref_level_3D(0,0); CHKERRQ(ierr);
129 
130  // stl::bitset see for more details
131  BitRefLevel bit_level0;
132  bit_level0.set(0);
133  EntityHandle meshset_level0;
134  rval = moab.create_meshset(MESHSET_SET,meshset_level0); CHKERRQ_MOAB(rval);
135  ierr = m_field.seed_ref_level_3D(0,bit_level0); CHKERRQ(ierr);
136  ierr = m_field.get_entities_by_ref_level(bit_level0,BitRefLevel().set(),meshset_level0); CHKERRQ(ierr);
137 
138  /***/
139  //Define problem
140 
141  //Fields
142  int field_rank=3;
143  ierr = m_field.add_field("DISPLACEMENT",H1,AINSWORTH_LEGENDRE_BASE,field_rank); CHKERRQ(ierr);
144  ierr = m_field.add_field("MESH_NODE_POSITIONS",H1,AINSWORTH_LEGENDRE_BASE,field_rank,MB_TAG_SPARSE,MF_ZERO); CHKERRQ(ierr);
145  ierr = m_field.add_field("LAGRANGE_MUL_TRAC",NOFIELD,NOBASE,6); CHKERRQ(ierr);
146  ierr = m_field.add_field("LAGRANGE_MUL_RIGID_TRANS",NOFIELD,NOBASE,3); CHKERRQ(ierr); //Control 3 rigid body translations in x, y and z axis
147  ierr = m_field.add_field("LAGRANGE_MUL_RIGID_ROT",NOFIELD,NOBASE,3); CHKERRQ(ierr); //Controla 3 rigid body rotations about x, y and z axis
148 
149 
150  //add entitities (by tets) to the field
151  ierr = m_field.add_ents_to_field_by_type(0,MBTET,"DISPLACEMENT"); CHKERRQ(ierr);
152  ierr = m_field.add_ents_to_field_by_type(0,MBTET,"MESH_NODE_POSITIONS"); CHKERRQ(ierr);
153 
154  //see Hierarchic Finite Element Bases on Unstructured Tetrahedral Meshes (Mark Ainsworth & Joe Coyle)
155  //int order = 5;
156  ierr = m_field.set_field_order(0,MBTET,"DISPLACEMENT",order); CHKERRQ(ierr);
157  ierr = m_field.set_field_order(0,MBTRI,"DISPLACEMENT",order); CHKERRQ(ierr);
158  ierr = m_field.set_field_order(0,MBEDGE,"DISPLACEMENT",order); CHKERRQ(ierr);
159  ierr = m_field.set_field_order(0,MBVERTEX,"DISPLACEMENT",1); CHKERRQ(ierr);
160 
161  ierr = m_field.set_field_order(0,MBTET,"MESH_NODE_POSITIONS",2); CHKERRQ(ierr);
162  ierr = m_field.set_field_order(0,MBTRI,"MESH_NODE_POSITIONS",2); CHKERRQ(ierr);
163  ierr = m_field.set_field_order(0,MBEDGE,"MESH_NODE_POSITIONS",2); CHKERRQ(ierr);
164  ierr = m_field.set_field_order(0,MBVERTEX,"MESH_NODE_POSITIONS",1); CHKERRQ(ierr);
165 
166 
167  //FE
168  //Define FE
169  //define eleatic element
170  boost::shared_ptr<Hooke<adouble> > hooke_adouble(new Hooke<adouble>);
171  boost::shared_ptr<Hooke<double> > hooke_double(new Hooke<double>);
172  NonlinearElasticElement elastic(m_field,2);
173  ierr = elastic.setBlocks(hooke_double,hooke_adouble); CHKERRQ(ierr);
174  ierr = elastic.addElement("ELASTIC","DISPLACEMENT"); CHKERRQ(ierr);
175  ierr = elastic.setOperators("DISPLACEMENT","MESH_NODE_POSITIONS",false,true); CHKERRQ(ierr);
176 
177 
178  BCs_RVELagrange_Trac lagrangian_element(m_field);
179  lagrangian_element.addLagrangiangElement("LAGRANGE_ELEM","DISPLACEMENT","LAGRANGE_MUL_TRAC","MESH_NODE_POSITIONS");
180  lagrangian_element.addLagrangiangElement("LAGRANGE_ELEM_TRANS","DISPLACEMENT","LAGRANGE_MUL_RIGID_TRANS","MESH_NODE_POSITIONS");
181  lagrangian_element.addLagrangiangElement("LAGRANGE_ELEM_ROT","DISPLACEMENT","LAGRANGE_MUL_RIGID_ROT","MESH_NODE_POSITIONS");
182 
183  //define problems
184  ierr = m_field.add_problem("ELASTIC_MECHANICS"); CHKERRQ(ierr);
185 
186  //set finite elements for problem
187  ierr = m_field.modify_problem_add_finite_element("ELASTIC_MECHANICS","ELASTIC"); CHKERRQ(ierr);
188  ierr = m_field.modify_problem_add_finite_element("ELASTIC_MECHANICS","LAGRANGE_ELEM"); CHKERRQ(ierr);
189  ierr = m_field.modify_problem_add_finite_element("ELASTIC_MECHANICS","LAGRANGE_ELEM_TRANS"); CHKERRQ(ierr);
190  ierr = m_field.modify_problem_add_finite_element("ELASTIC_MECHANICS","LAGRANGE_ELEM_ROT"); CHKERRQ(ierr);
191 
192 
193  //set refinement level for problem
194  ierr = m_field.modify_problem_ref_level_add_bit("ELASTIC_MECHANICS",bit_level0); CHKERRQ(ierr);
195 
196  /***/
197  //Declare problem
198  /****/
199  //build database
200 
201  //build field
202  ierr = m_field.build_fields(); CHKERRQ(ierr);
203  Projection10NodeCoordsOnField ent_method_material(m_field,"MESH_NODE_POSITIONS");
204  ierr = m_field.loop_dofs("MESH_NODE_POSITIONS",ent_method_material); CHKERRQ(ierr);
205 
206  //build finite elemnts
207  ierr = m_field.build_finite_elements(); CHKERRQ(ierr);
208 
209  //build adjacencies
210  ierr = m_field.build_adjacencies(bit_level0); CHKERRQ(ierr);
211 
212  //build problem
213  ierr = m_field.build_problems(); CHKERRQ(ierr);
214 
215  /****/
216  //mesh partitioning
217 
218  //partition
219  ierr = m_field.partition_problem("ELASTIC_MECHANICS"); CHKERRQ(ierr);
220  ierr = m_field.partition_finite_elements("ELASTIC_MECHANICS"); CHKERRQ(ierr);
221  //what are ghost nodes, see Petsc Manual
222  ierr = m_field.partition_ghost_dofs("ELASTIC_MECHANICS"); CHKERRQ(ierr);
223 
224  //print bcs
225  ierr = m_field.print_cubit_displacement_set(); CHKERRQ(ierr);
226  ierr = m_field.print_cubit_force_set(); CHKERRQ(ierr);
227  //print block sets with materials
228  ierr = m_field.print_cubit_materials_set(); CHKERRQ(ierr);
229 
230  //create matrices (here F, D and Aij are matrices for the full problem)
231  vector<Vec> F(6);
232  ierr = m_field.VecCreateGhost("ELASTIC_MECHANICS",ROW,&F[0]); CHKERRQ(ierr);
233  for(int ii = 1;ii<6;ii++) {
234  ierr = VecDuplicate(F[0],&F[ii]); CHKERRQ(ierr);
235  }
236 
237  Vec D;
238  ierr = m_field.VecCreateGhost("ELASTIC_MECHANICS",COL,&D); CHKERRQ(ierr);
239 
240  Mat Aij;
241  ierr = m_field.MatCreateMPIAIJWithArrays("ELASTIC_MECHANICS",&Aij); CHKERRQ(ierr);
242 
243  ierr = VecZeroEntries(D); CHKERRQ(ierr);
244  ierr = VecGhostUpdateBegin(D,INSERT_VALUES,SCATTER_FORWARD); CHKERRQ(ierr);
245  ierr = VecGhostUpdateEnd(D,INSERT_VALUES,SCATTER_FORWARD); CHKERRQ(ierr);
246  ierr = m_field.set_global_ghost_vector(
247  "ELASTIC_MECHANICS",ROW,D,INSERT_VALUES,SCATTER_REVERSE
248  ); CHKERRQ(ierr);
249 
250  for(int ii = 0;ii<6;ii++) {
251  ierr = VecZeroEntries(F[ii]); CHKERRQ(ierr);
252  ierr = VecGhostUpdateBegin(F[ii],INSERT_VALUES,SCATTER_FORWARD); CHKERRQ(ierr);
253  ierr = VecGhostUpdateEnd(F[ii],INSERT_VALUES,SCATTER_FORWARD); CHKERRQ(ierr);
254  }
255  ierr = MatZeroEntries(Aij); CHKERRQ(ierr);
256 
257  //Assemble F and Aij
258  elastic.getLoopFeLhs().snes_B = Aij;
259  ierr = m_field.loop_finite_elements("ELASTIC_MECHANICS","ELASTIC",elastic.getLoopFeLhs()); CHKERRQ(ierr);
260  lagrangian_element.setRVEBCsOperators("DISPLACEMENT","LAGRANGE_MUL_TRAC","MESH_NODE_POSITIONS",Aij,F);
261  ierr = m_field.loop_finite_elements("ELASTIC_MECHANICS","LAGRANGE_ELEM",lagrangian_element.getLoopFeRVEBCLhs()); CHKERRQ(ierr);
262  ierr = m_field.loop_finite_elements("ELASTIC_MECHANICS","LAGRANGE_ELEM",lagrangian_element.getLoopFeRVEBCRhs()); CHKERRQ(ierr);
263 
264  BCs_RVELagrange_Trac_Rigid_Trans lagrangian_element_rigid_body_trans(m_field);
265  lagrangian_element_rigid_body_trans.setRVEBCsRigidBodyTranOperators("DISPLACEMENT","LAGRANGE_MUL_RIGID_TRANS",Aij, lagrangian_element.setOfRVEBC);
266  ierr = m_field.loop_finite_elements("ELASTIC_MECHANICS","LAGRANGE_ELEM_TRANS",lagrangian_element_rigid_body_trans.getLoopFeRVEBCLhs()); CHKERRQ(ierr);
267 
268  BCs_RVELagrange_Trac_Rigid_Rot lagrangian_element_rigid_body_rot(m_field);
269  lagrangian_element_rigid_body_rot.setRVEBCsRigidBodyRotOperators("DISPLACEMENT","LAGRANGE_MUL_RIGID_ROT",Aij, lagrangian_element.setOfRVEBC);
270  ierr = m_field.loop_finite_elements("ELASTIC_MECHANICS","LAGRANGE_ELEM_ROT",lagrangian_element_rigid_body_rot.getLoopFeRVEBCLhs()); CHKERRQ(ierr);
271 
272 // //Matrix View
273 // MatView(Aij,PETSC_VIEWER_DRAW_WORLD);//PETSC_VIEWER_STDOUT_WORLD);
274 // std::string wait;
275 // std::cin >> wait;
276 
277  ierr = MatAssemblyBegin(Aij,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
278  ierr = MatAssemblyEnd(Aij,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
279 
280  for(int ii = 0;ii<6;ii++) {
281  ierr = VecGhostUpdateBegin(F[ii],ADD_VALUES,SCATTER_REVERSE); CHKERRQ(ierr);
282  ierr = VecGhostUpdateEnd(F[ii],ADD_VALUES,SCATTER_REVERSE); CHKERRQ(ierr);
283  ierr = VecAssemblyBegin(F[ii]); CHKERRQ(ierr);
284  ierr = VecAssemblyEnd(F[ii]); CHKERRQ(ierr);
285  }
286 
287 // //ierr = VecView(F,PETSC_VIEWER_STDOUT_WORLD); CHKERRQ(ierr);
288 // //ierr = MatView(Aij,PETSC_VIEWER_STDOUT_WORLD); CHKERRQ(ierr);
289 //
290 ////Matrix View
291 //MatView(Aij,PETSC_VIEWER_DRAW_WORLD);//PETSC_VIEWER_STDOUT_WORLD);
292 //std::string wait;
293 //std::cin >> wait;
294 
295 
296 
297  // Volume calculation
298  //==============================================================================================================================
299  VolumeElementForcesAndSourcesCore vol_elem(m_field);
300  Vec volume_vec;
301  int volume_vec_ghost[] = { 0 };
302  ierr = VecCreateGhost(
303  PETSC_COMM_WORLD,(!m_field.get_comm_rank())?1:0,1,1,volume_vec_ghost,&volume_vec
304  ); CHKERRQ(ierr);
305  ierr = VecZeroEntries(volume_vec); CHKERRQ(ierr);
306  vol_elem.getOpPtrVector().push_back(new VolumeCalculation("DISPLACEMENT",volume_vec));
307 
308  ierr = m_field.loop_finite_elements("ELASTIC_MECHANICS","ELASTIC", vol_elem); CHKERRQ(ierr);
309 
310  ierr = VecAssemblyBegin(volume_vec); CHKERRQ(ierr);
311  ierr = VecAssemblyEnd(volume_vec); CHKERRQ(ierr);
312  double rve_volume;
313  ierr = VecSum(volume_vec,&rve_volume); CHKERRQ(ierr);
314 
315  ierr = VecView(volume_vec,PETSC_VIEWER_STDOUT_WORLD); CHKERRQ(ierr);
316  ierr = PetscPrintf(PETSC_COMM_WORLD,"RVE Volume %3.2g\n",rve_volume); CHKERRQ(ierr);
317  ierr = VecDestroy(&volume_vec);
318 
319  //==============================================================================================================================
320  //Post processing
321  //==============================================================================================================================
322  struct MyPostProc: public PostProcVolumeOnRefinedMesh {
323  bool doPreProcess;
324  bool doPostProcess;
325  MyPostProc(MoFEM::Interface &m_field):
326  PostProcVolumeOnRefinedMesh(m_field),
327  doPreProcess(true),
328  doPostProcess(true)
329  {}
330  void setDoPreProcess() { doPreProcess = true; }
331  void unSetDoPreProcess() { doPreProcess = false; }
332  void setDoPostProcess() { doPostProcess = true; }
333  void unSetDoPostProcess() { doPostProcess = false; }
334 
335 
336 
337  PetscErrorCode preProcess() {
338  PetscFunctionBegin;
339  if(doPreProcess) {
340  ierr = PostProcVolumeOnRefinedMesh::preProcess(); CHKERRQ(ierr);
341  }
342  PetscFunctionReturn(0);
343  }
344  PetscErrorCode postProcess() {
345  PetscFunctionBegin;
346  if(doPostProcess) {
347  ierr = PostProcVolumeOnRefinedMesh::postProcess(); CHKERRQ(ierr);
348  }
349  PetscFunctionReturn(0);
350  }
351  };
352 
353  MyPostProc post_proc(m_field);
354 
355  ierr = post_proc.generateReferenceElementMesh(); CHKERRQ(ierr);
356  ierr = post_proc.addFieldValuesPostProc("DISPLACEMENT"); CHKERRQ(ierr);
357  ierr = post_proc.addFieldValuesGradientPostProc("DISPLACEMENT"); CHKERRQ(ierr);
358  ierr = post_proc.addFieldValuesPostProc("MESH_NODE_POSITIONS"); CHKERRQ(ierr);
359  ierr = post_proc.addFieldValuesGradientPostProc("MESH_NODE_POSITIONS"); CHKERRQ(ierr);
360 
361  for(
362  map<int,NonlinearElasticElement::BlockData>::iterator sit = elastic.setOfBlocks.begin();
363  sit != elastic.setOfBlocks.end(); sit++
364  ) {
365  post_proc.getOpPtrVector().push_back(
366  new PostPorcStress(
367  post_proc.postProcMesh,
368  post_proc.mapGaussPts,
369  "DISPLACEMENT",
370  sit->second,
371  post_proc.commonData,
372  false
373  )
374  );
375  }
376  ierr = VecGhostUpdateBegin(D,INSERT_VALUES,SCATTER_FORWARD); CHKERRQ(ierr);
377  ierr = VecGhostUpdateEnd(D,INSERT_VALUES,SCATTER_FORWARD); CHKERRQ(ierr);
378  ierr = m_field.set_global_ghost_vector(
379  "ELASTIC_MECHANICS",ROW,D,INSERT_VALUES,SCATTER_REVERSE
380  ); CHKERRQ(ierr);
381 
382  //==============================================================================================================================
383 
384  //Solver
385  KSP solver;
386  ierr = KSPCreate(PETSC_COMM_WORLD,&solver); CHKERRQ(ierr);
387  ierr = KSPSetOperators(solver,Aij,Aij); CHKERRQ(ierr);
388  ierr = KSPSetFromOptions(solver); CHKERRQ(ierr);
389  ierr = KSPSetUp(solver); CHKERRQ(ierr);
390 
391  MatrixDouble Dmat;
392  Dmat.resize(6,6);
393  Dmat.clear();
394 
395  //calculate honmogenised stress
396  //create a vector for 6 components of homogenized stress
397  Vec stress_homo;
398  int stress_homo_ghost[] = { 0,1,2,3,4,5,6 };
399  ierr = VecCreateGhost(
400  PETSC_COMM_WORLD,(!m_field.get_comm_rank())?6:0,6,6,stress_homo_ghost,&stress_homo
401  ); CHKERRQ(ierr);
402 
403  lagrangian_element.setRVEBCsHomoStressOperators("DISPLACEMENT","LAGRANGE_MUL_TRAC","MESH_NODE_POSITIONS",stress_homo);
404 
405  PetscScalar *avec;
406  ierr = VecGetArray(stress_homo,&avec); CHKERRQ(ierr);
407  for(int ii = 0;ii<6;ii++) {
408  ierr = VecZeroEntries(D); CHKERRQ(ierr);
409  ierr = KSPSolve(solver,F[ii],D); CHKERRQ(ierr);
410  ierr = VecGhostUpdateBegin(D,INSERT_VALUES,SCATTER_FORWARD); CHKERRQ(ierr);
411  ierr = VecGhostUpdateEnd(D,INSERT_VALUES,SCATTER_FORWARD); CHKERRQ(ierr);
412  ierr = m_field.set_global_ghost_vector(
413  "ELASTIC_MECHANICS",ROW,D,INSERT_VALUES,SCATTER_REVERSE
414  ); CHKERRQ(ierr);
415 
416  //post-processing the resutls
417  post_proc.setDoPreProcess();
418  post_proc.unSetDoPostProcess();
419  ierr = m_field.loop_finite_elements(
420  "ELASTIC_MECHANICS","ELASTIC",post_proc
421  ); CHKERRQ(ierr);
422  post_proc.unSetDoPreProcess();
423  post_proc.setDoPostProcess();
424  {
425  ostringstream sss;
426  sss << "mode_" << "Disp" << "_" << ii << ".h5m";
427  rval = post_proc.postProcMesh.write_file(
428  sss.str().c_str(),"MOAB","PARALLEL=WRITE_PART"
429  ); CHKERRQ_MOAB(rval);
430  }
431 
432  ierr = VecZeroEntries(stress_homo); CHKERRQ(ierr);
433  ierr = m_field.loop_finite_elements(
434  "ELASTIC_MECHANICS","LAGRANGE_ELEM",lagrangian_element.getLoopFeRVEBCStress()
435  ); CHKERRQ(ierr);
436  ierr = PetscOptionsGetReal(
437  PETSC_NULL,"-my_rve_volume",&rve_volume,PETSC_NULL
438  ); CHKERRQ(ierr);
439  ierr = VecAssemblyBegin(stress_homo); CHKERRQ(ierr);
440  ierr = VecAssemblyEnd(stress_homo); CHKERRQ(ierr);
441  ierr = VecScale(stress_homo,1.0/rve_volume); CHKERRQ(ierr);
442  ierr = VecGhostUpdateBegin(stress_homo,INSERT_VALUES,SCATTER_FORWARD); CHKERRQ(ierr);
443  ierr = VecGhostUpdateEnd(stress_homo,INSERT_VALUES,SCATTER_FORWARD); CHKERRQ(ierr);
444  for(int jj=0; jj<6; jj++) {
445  Dmat(jj,ii) = avec[jj];
446  }
447  }
448 
449  ierr = VecRestoreArray(stress_homo,&avec); CHKERRQ(ierr);
450  PetscPrintf(
451  PETSC_COMM_WORLD,"\nHomogenised Stiffens Matrix = \n\n"
452  );
453 
454  for(int ii=0; ii<6; ii++) {
455  PetscPrintf(
456  PETSC_COMM_WORLD,
457  "stress %d\t\t%8.5e\t\t%8.5e\t\t%8.5e\t\t%8.5e\t\t%8.5e\t\t%8.5e\n",
458  ii,Dmat(ii,0),Dmat(ii,1),Dmat(ii,2),Dmat(ii,3),Dmat(ii,4),Dmat(ii,5)
459  );
460  }
461  //==============================================================================================================================
462 
463  //Open mesh_file_name.txt for writing
464  ofstream myfile;
465  myfile.open ((string(mesh_file_name)+".txt").c_str());
466 
467  //Output displacements
468  myfile << "<<<< Homonenised stress >>>>>" << endl;
469 
470  if(pcomm->rank()==0){
471  for(int ii=0; ii<6; ii++){
472  myfile << boost::format("%.3lf") % roundn(Dmat(ii,0)) << endl;
473  avec++;
474  }
475  }
476  //Close mesh_file_name.txt
477  myfile.close();
478 
479 
480 
481  //detroy matrices
482  for(int ii = 0;ii<6;ii++) {
483  ierr = VecDestroy(&F[ii]); CHKERRQ(ierr);
484  }
485  ierr = VecDestroy(&D); CHKERRQ(ierr);
486  ierr = MatDestroy(&Aij); CHKERRQ(ierr);
487  ierr = KSPDestroy(&solver); CHKERRQ(ierr);
488  ierr = VecDestroy(&stress_homo); CHKERRQ(ierr);
489 
490  ierr = PetscTime(&v2);CHKERRQ(ierr);
491  ierr = PetscGetCPUTime(&t2);CHKERRQ(ierr);
492 
493  PetscSynchronizedPrintf(PETSC_COMM_WORLD,"Total Rank %d Time = %f CPU Time = %f\n",pcomm->rank(),v2-v1,t2-t1);
494 
495  }
496  CATCH_ERRORS;
497 
498  MoFEM::Core::Finalize();
499 
500 }
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
BCs_RVELagrange_Trac_Rigid_Rot
Definition: BCs_RVELagrange_Trac_Rigid_Rot.hpp:22
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.
MoFEM::CoreTmp< 0 >
Core (interface) class.
Definition: Core.hpp:82
H1
@ H1
continuous field
Definition: definitions.h:85
MyPostProc
Definition: plot_base.cpp:34
BCs_RVELagrange_Disp.hpp
EntityHandle
CHKERRQ_MOAB
#define CHKERRQ_MOAB(a)
check error code of MoAB function
Definition: definitions.h:467
NOBASE
@ NOBASE
Definition: definitions.h:59
NonlinearElasticElement
structure grouping operators and data used for calculation of nonlinear elastic element
Definition: HookeElement.hpp:27
PostProcTemplateVolumeOnRefinedMesh< MoFEM::VolumeElementForcesAndSourcesCore >::postProcess
MoFEMErrorCode postProcess()
function is run at the end of loop
Definition: PostProcOnRefMesh.hpp:797
MoFEM::Types::MatrixDouble
UBlasMatrix< double > MatrixDouble
Definition: Types.hpp:77
MoFEM::CoreInterface::get_comm_rank
virtual int get_comm_rank() const =0
BCs_RVELagrange_Disp::getLoopFeRVEBCRhs
MyTriFE & getLoopFeRVEBCRhs()
Definition: BCs_RVELagrange_Disp.hpp:36
MoFEM.hpp
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
RND_EPS
#define RND_EPS
Definition: homogenisation_trac_atom_test.cpp:52
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.
order
constexpr int order
Definition: dg_projection.cpp:18
MoFEM::DeprecatedCoreInterface
Deprecated interface functions.
Definition: DeprecatedCoreInterface.hpp:16
ROW
@ ROW
Definition: definitions.h:136
MoFEM::Interface
DeprecatedCoreInterface Interface
Definition: Interface.hpp:2010
MoFEM::Exceptions::rval
static MoFEMErrorCodeGeneric< moab::ErrorCode > rval
Definition: Exceptions.hpp:74
PostProcStresses.hpp
Post-processing stresses for non-linear analysis.
PostProcTemplateVolumeOnRefinedMesh< MoFEM::VolumeElementForcesAndSourcesCore >::preProcess
MoFEMErrorCode preProcess()
function is run at the beginning of loop
Definition: PostProcOnRefMesh.hpp:696
MoFEM
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:10
MoFEM::PetscOptionsGetReal
PetscErrorCode PetscOptionsGetReal(PetscOptions *, const char pre[], const char name[], PetscReal *dval, PetscBool *set)
Definition: DeprecatedPetsc.hpp:152
BCs_RVELagrange_Disp::getLoopFeRVEBCLhs
MyTriFE & getLoopFeRVEBCLhs()
Definition: BCs_RVELagrange_Disp.hpp:35
MoFEM::CoreInterface::build_finite_elements
virtual MoFEMErrorCode build_finite_elements(int verb=DEFAULT_VERBOSITY)=0
Build finite elements.
Hooke.hpp
Implementation of linear elastic material.
BCs_RVELagrange_Trac_Rigid_Trans.hpp
COL
@ COL
Definition: definitions.h:136
NonLinearElasticElement.hpp
Operators and data structures for non-linear elastic analysis.
Projection10NodeCoordsOnField.hpp
BCs_RVELagrange_Trac
Definition: BCs_RVELagrange_Trac.hpp:22
PostProcOnRefMesh.hpp
Post-process fields on refined mesh.
NonlinearElasticElement::addElement
MoFEMErrorCode addElement(const std::string element_name, const std::string spatial_position_field_name, const std::string material_position_field_name="MESH_NODE_POSITIONS", const bool ale=false)
Definition: NonLinearElasticElement.cpp:1120
VolumeCalculation.hpp
Operator can be used with any volume element to calculate sum of volumes of all volumes in the set.
MyPostProc::generateReferenceElementMesh
MoFEMErrorCode generateReferenceElementMesh()
Definition: plot_base.cpp:427
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
MoFEM::VolumeElementForcesAndSourcesCore
Volume finite element base.
Definition: VolumeElementForcesAndSourcesCore.hpp:26
main
int main(int argc, char *argv[])
Definition: homogenisation_trac_atom_test.cpp:86
MoFEM::CoreInterface::modify_problem_ref_level_add_bit
virtual MoFEMErrorCode modify_problem_ref_level_add_bit(const std::string &name_problem, const BitRefLevel &bit)=0
add ref level to problem
NonlinearElasticElement::setOperators
MoFEMErrorCode setOperators(const std::string spatial_position_field_name, const std::string material_position_field_name="MESH_NODE_POSITIONS", const bool ale=false, const bool field_disp=false)
Set operators to calculate left hand tangent matrix and right hand residual.
Definition: NonLinearElasticElement.cpp:1153
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:22
BCs_RVELagrange_Trac::setRVEBCsOperators
PetscErrorCode setRVEBCsOperators(string field_name, string lagrang_field_name, string mesh_nodals_positions, Mat aij, vector< Vec > &f)
Definition: BCs_RVELagrange_Trac.hpp:383
convert.n
n
Definition: convert.py:82
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
CATCH_ERRORS
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:385
PostProcVolumeOnRefinedMesh
Post processing.
Definition: PostProcOnRefMesh.hpp:955
PostPorcStress
DEPRECATED typedef PostProcStress PostPorcStress
Definition: PostProcStresses.hpp:193
BCs_RVELagrange_Trac_Rigid_Rot.hpp
MoFEM::Core
CoreTmp< 0 > Core
Definition: Core.hpp:1148
help
static char help[]
Definition: homogenisation_trac_atom_test.cpp:51
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
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
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
BCs_RVELagrange_Trac_Rigid_Rot::setRVEBCsRigidBodyRotOperators
PetscErrorCode setRVEBCsRigidBodyRotOperators(string field_name, string lagrang_field_name, Mat aij, map< int, RVEBC_Data > setOfRVEBC)
Definition: BCs_RVELagrange_Trac_Rigid_Rot.hpp:105
MoFEM::CoreInterface::modify_problem_add_finite_element
virtual MoFEMErrorCode modify_problem_add_finite_element(const std::string name_problem, const std::string &fe_name)=0
add finite element to problem, this add entities assigned to finite element to a particular problem
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::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
BCs_RVELagrange_Trac::setRVEBCsHomoStressOperators
PetscErrorCode setRVEBCsHomoStressOperators(string field_name, string lagrang_field_name, string mesh_nodals_positions, Vec stress_homo)
Definition: BCs_RVELagrange_Trac.hpp:448
BCs_RVELagrange_Trac.hpp
MoFEM::CoreInterface::add_problem
virtual MoFEMErrorCode add_problem(const std::string &name, enum MoFEMTypes bh=MF_EXCL, int verb=DEFAULT_VERBOSITY)=0
Add problem.
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.
NonlinearElasticElement::setBlocks
MoFEMErrorCode setBlocks(boost::shared_ptr< FunctionsToCalculatePiolaKirchhoffI< double >> materialDoublePtr, boost::shared_ptr< FunctionsToCalculatePiolaKirchhoffI< adouble >> materialAdoublePtr)
Definition: NonLinearElasticElement.cpp:1086
F
@ F
Definition: free_surface.cpp:394
Hooke
Hook equation.
Definition: Hooke.hpp:19
NOFIELD
@ NOFIELD
scalar or vector of scalars describe (no true field)
Definition: definitions.h:84
roundn
double roundn(double n)
Definition: homogenisation_trac_atom_test.cpp:55
Generated by Doxygen 1.8.17 and hosted at University of Glasgow