v0.14.0
homogenisation_disp_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 <BCs_RVELagrange_Disp.hpp>
38 #include <VolumeCalculation.hpp>
39 
40 #include <boost/ptr_container/ptr_map.hpp>
41 #include <PostProcOnRefMesh.hpp>
42 #include <PostProcStresses.hpp>
43 
44 
45 
46 
47 static char help[] = "...\n\n";
48 #define RND_EPS 1e-6
49 
50 
51 //Rounding
52 double roundn(double n)
53 {
54  //break n into fractional part (fract) and integral part (intp)
55  double fract, intp;
56  fract = modf(n,&intp);
57 
58  // //round up
59  // if (fract>=.5)
60  // {
61  // n*=10;
62  // ceil(n);
63  // n/=10;
64  // }
65  // //round down
66  // if (fract<=.5)
67  // {
68  // n*=10;
69  // floor(n);
70  // n/=10;
71  // }
72  // case where n approximates zero, set n to "positive" zero
73  if (abs(intp)==0)
74  {
75  if(abs(fract)<=RND_EPS)
76  {
77  n=0.000;
78  }
79  }
80  return n;
81 }
82 
83 
84 int main(int argc, char *argv[]) {
85 
86  MoFEM::Core::Initialize(&argc,&argv,(char *)0,help);
87 
88  try {
89 
90  moab::Core mb_instance;
91  moab::Interface& moab = mb_instance;
92  int rank;
93  MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
94 
95  //Reade parameters from line command
96  PetscBool flg = PETSC_TRUE;
97  char mesh_file_name[255];
98  ierr = PetscOptionsGetString(PETSC_NULL,"-my_file",mesh_file_name,255,&flg); CHKERRQ(ierr);
99  if(flg != PETSC_TRUE) {
100  SETERRQ(PETSC_COMM_SELF,1,"*** ERROR -my_file (MESH FILE NEEDED)");
101  }
102  PetscInt order;
103  ierr = PetscOptionsGetInt(PETSC_NULL,"-my_order",&order,&flg); CHKERRQ(ierr);
104  if(flg != PETSC_TRUE) {
105  order = 1;
106  }
107 
108  //Read mesh to MOAB
109  const char *option;
110  option = "";//"PARALLEL=BCAST;";//;DEBUG_IO";
111  rval = moab.load_file(mesh_file_name, 0, option); CHKERRQ_MOAB(rval);
112  ParallelComm* pcomm = ParallelComm::get_pcomm(&moab,MYPCOMM_INDEX);
113  if(pcomm == NULL) pcomm = new ParallelComm(&moab,PETSC_COMM_WORLD);
114 
115  //We need that for code profiling
116  PetscLogDouble t1,t2;
117  PetscLogDouble v1,v2;
118  ierr = PetscTime(&v1); CHKERRQ(ierr);
119  ierr = PetscGetCPUTime(&t1); CHKERRQ(ierr);
120 
121  //Create MoFEM (Joseph) database
122  MoFEM::Core core(moab);
123  MoFEM::Interface& m_field = core;
124 
125  //ref meshset ref level 0
126  ierr = m_field.seed_ref_level_3D(0,0); CHKERRQ(ierr);
127 
128  // stl::bitset see for more details
129  BitRefLevel bit_level0;
130  bit_level0.set(0);
131  EntityHandle meshset_level0;
132  rval = moab.create_meshset(MESHSET_SET,meshset_level0); CHKERRQ_MOAB(rval);
133  ierr = m_field.seed_ref_level_3D(0,bit_level0); CHKERRQ(ierr);
134  ierr = m_field.get_entities_by_ref_level(bit_level0,BitRefLevel().set(),meshset_level0); CHKERRQ(ierr);
135 
136  /***/
137  //Define problem
138 
139  //Fields
140  int field_rank=3;
141  ierr = m_field.add_field("DISPLACEMENT",H1,AINSWORTH_LEGENDRE_BASE,field_rank,MB_TAG_SPARSE,MF_ZERO); CHKERRQ(ierr);
142  ierr = m_field.add_field("MESH_NODE_POSITIONS",H1,AINSWORTH_LEGENDRE_BASE,field_rank,MB_TAG_SPARSE,MF_ZERO); CHKERRQ(ierr);
143  ierr = m_field.add_field("LAGRANGE_MUL_DISP",H1,AINSWORTH_LEGENDRE_BASE,field_rank); CHKERRQ(ierr);
144 
145  //Add entitities to the fields
146  //============================================================================================================
147  ierr = m_field.add_ents_to_field_by_type(0,MBTET,"DISPLACEMENT"); CHKERRQ(ierr);
148  ierr = m_field.add_ents_to_field_by_type(0,MBTET,"MESH_NODE_POSITIONS"); CHKERRQ(ierr);
149 
150  Range SurfacesFaces;
151  ierr = m_field.get_cubit_msId_entities_by_dimension(103,SIDESET,2,SurfacesFaces,true); CHKERRQ(ierr);
152  ierr = PetscPrintf(PETSC_COMM_WORLD,"number of SideSet 103 = %d\n",SurfacesFaces.size()); CHKERRQ(ierr);
153 
154  EntityHandle BoundFacesMeshset;
155  rval = moab.create_meshset(MESHSET_SET,BoundFacesMeshset); CHKERRQ_MOAB(rval);
156  rval = moab.add_entities(BoundFacesMeshset,SurfacesFaces); CHKERRQ_MOAB(rval);
157  ierr = m_field.seed_ref_level_MESHSET(BoundFacesMeshset,BitRefLevel().set()); CHKERRQ(ierr);
158  ierr = m_field.add_ents_to_field_by_type(BoundFacesMeshset,MBTRI,"LAGRANGE_MUL_DISP"); CHKERRQ(ierr);
159  //============================================================================================================
160 
161  //set approximaiton order
162  //see Hierarchic Finite Element Bases on Unstructured Tetrahedral Meshes (Mark Ainsworth & Joe Coyle)
163  ierr = m_field.set_field_order(0,MBTET,"DISPLACEMENT",order); CHKERRQ(ierr);
164  ierr = m_field.set_field_order(0,MBTRI,"DISPLACEMENT",order); CHKERRQ(ierr);
165  ierr = m_field.set_field_order(0,MBEDGE,"DISPLACEMENT",order); CHKERRQ(ierr);
166  ierr = m_field.set_field_order(0,MBVERTEX,"DISPLACEMENT",1); CHKERRQ(ierr);
167 
168  ierr = m_field.set_field_order(0,MBTET,"MESH_NODE_POSITIONS",2); CHKERRQ(ierr);
169  ierr = m_field.set_field_order(0,MBTRI,"MESH_NODE_POSITIONS",2); CHKERRQ(ierr);
170  ierr = m_field.set_field_order(0,MBEDGE,"MESH_NODE_POSITIONS",2); CHKERRQ(ierr);
171  ierr = m_field.set_field_order(0,MBVERTEX,"MESH_NODE_POSITIONS",1); CHKERRQ(ierr);
172 
173  ierr = m_field.set_field_order(0,MBTRI,"LAGRANGE_MUL_DISP",order); CHKERRQ(ierr);
174  ierr = m_field.set_field_order(0,MBEDGE,"LAGRANGE_MUL_DISP",order); CHKERRQ(ierr);
175  ierr = m_field.set_field_order(0,MBVERTEX,"LAGRANGE_MUL_DISP",1); CHKERRQ(ierr);
176 
177 
178  //Define FE
179  //define eleatic element
180  boost::shared_ptr<Hooke<adouble> > hooke_adouble(new Hooke<adouble>);
181  boost::shared_ptr<Hooke<double> > hooke_double(new Hooke<double>);
182  NonlinearElasticElement elastic(m_field,2);
183  ierr = elastic.setBlocks(hooke_double,hooke_adouble); CHKERRQ(ierr);
184  ierr = elastic.addElement("ELASTIC","DISPLACEMENT"); CHKERRQ(ierr);
185  ierr = elastic.setOperators("DISPLACEMENT","MESH_NODE_POSITIONS",false,true); CHKERRQ(ierr);
186 
187 
188  BCs_RVELagrange_Disp lagrangian_element(m_field);
189  lagrangian_element.addLagrangiangElement("LAGRANGE_ELEM","DISPLACEMENT","LAGRANGE_MUL_DISP","MESH_NODE_POSITIONS");
190 
191  //define problems
192  ierr = m_field.add_problem("ELASTIC_MECHANICS"); CHKERRQ(ierr);
193 
194  //set finite elements for problem
195  ierr = m_field.modify_problem_add_finite_element("ELASTIC_MECHANICS","ELASTIC"); CHKERRQ(ierr);
196  ierr = m_field.modify_problem_add_finite_element("ELASTIC_MECHANICS","LAGRANGE_ELEM"); CHKERRQ(ierr);
197 
198 
199  //set refinement level for problem
200  ierr = m_field.modify_problem_ref_level_add_bit("ELASTIC_MECHANICS",bit_level0); CHKERRQ(ierr);
201 
202  /***/
203  //Declare problem
204 
205  /****/
206  //build database
207 
208  //build field
209  ierr = m_field.build_fields(); CHKERRQ(ierr);
210  Projection10NodeCoordsOnField ent_method_material(m_field,"MESH_NODE_POSITIONS");
211  ierr = m_field.loop_dofs("MESH_NODE_POSITIONS",ent_method_material); CHKERRQ(ierr);
212 
213  //build finite elemnts
214  ierr = m_field.build_finite_elements(); CHKERRQ(ierr);
215 
216  //build adjacencies
217  ierr = m_field.build_adjacencies(bit_level0); CHKERRQ(ierr);
218 
219  //build problem
220  ierr = m_field.build_problems(); CHKERRQ(ierr);
221 
222 
223  /****/
224  //mesh partitioning
225 
226  //partition
227  ierr = m_field.partition_problem("ELASTIC_MECHANICS"); CHKERRQ(ierr);
228  ierr = m_field.partition_finite_elements("ELASTIC_MECHANICS"); CHKERRQ(ierr);
229  //what are ghost nodes, see Petsc Manual
230  ierr = m_field.partition_ghost_dofs("ELASTIC_MECHANICS"); CHKERRQ(ierr);
231 
232  //print bcs
233  ierr = m_field.print_cubit_displacement_set(); CHKERRQ(ierr);
234  ierr = m_field.print_cubit_force_set(); CHKERRQ(ierr);
235  //print block sets with materials
236  ierr = m_field.print_cubit_materials_set(); CHKERRQ(ierr);
237 
238  //create matrices (here F, D and Aij are matrices for the full problem)
239  vector<Vec> F(6);
240  ierr = m_field.VecCreateGhost("ELASTIC_MECHANICS",ROW,&F[0]); CHKERRQ(ierr);
241  for(int ii = 1;ii<6;ii++) {
242  ierr = VecDuplicate(F[0],&F[ii]); CHKERRQ(ierr);
243  }
244 
245  Vec D;
246  ierr = m_field.VecCreateGhost("ELASTIC_MECHANICS",COL,&D); CHKERRQ(ierr);
247 
248  Mat Aij;
249  ierr = m_field.MatCreateMPIAIJWithArrays("ELASTIC_MECHANICS",&Aij); CHKERRQ(ierr);
250 
251  ierr = VecZeroEntries(D); CHKERRQ(ierr);
252  ierr = VecGhostUpdateBegin(D,INSERT_VALUES,SCATTER_FORWARD); CHKERRQ(ierr);
253  ierr = VecGhostUpdateEnd(D,INSERT_VALUES,SCATTER_FORWARD); CHKERRQ(ierr);
254  ierr = m_field.set_global_ghost_vector(
255  "ELASTIC_MECHANICS",ROW,D,INSERT_VALUES,SCATTER_REVERSE
256  ); CHKERRQ(ierr);
257 
258  for(int ii = 0;ii<6;ii++) {
259  ierr = VecZeroEntries(F[ii]); CHKERRQ(ierr);
260  ierr = VecGhostUpdateBegin(F[ii],INSERT_VALUES,SCATTER_FORWARD); CHKERRQ(ierr);
261  ierr = VecGhostUpdateEnd(F[ii],INSERT_VALUES,SCATTER_FORWARD); CHKERRQ(ierr);
262  }
263  ierr = MatZeroEntries(Aij); CHKERRQ(ierr);
264 
265  //Assemble F and Aij
266  elastic.getLoopFeLhs().snes_B = Aij;
267  ierr = m_field.loop_finite_elements("ELASTIC_MECHANICS","ELASTIC",elastic.getLoopFeLhs()); CHKERRQ(ierr);
268 
269  lagrangian_element.setRVEBCsOperators("DISPLACEMENT","LAGRANGE_MUL_DISP","MESH_NODE_POSITIONS",Aij,F);
270  ierr = m_field.loop_finite_elements("ELASTIC_MECHANICS","LAGRANGE_ELEM",lagrangian_element.getLoopFeRVEBCLhs()); CHKERRQ(ierr);
271  cout<<"after lagrangian_element.getLoopFeRVEBCLhs "<<endl;
272  ierr = m_field.loop_finite_elements("ELASTIC_MECHANICS","LAGRANGE_ELEM",lagrangian_element.getLoopFeRVEBCRhs()); CHKERRQ(ierr);
273  cout<<"after lagrangian_element.getLoopFeRVEBCRhs "<<endl;
274 //
275 // //Matrix View
276 // MatView(Aij,PETSC_VIEWER_DRAW_WORLD);//PETSC_VIEWER_STDOUT_WORLD);
277 // std::string wait;
278 // std::cin >> wait;
279 
280  ierr = MatAssemblyBegin(Aij,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
281  ierr = MatAssemblyEnd(Aij,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
282 
283  for(int ii = 0;ii<6;ii++) {
284  ierr = VecGhostUpdateBegin(F[ii],ADD_VALUES,SCATTER_REVERSE); CHKERRQ(ierr);
285  ierr = VecGhostUpdateEnd(F[ii],ADD_VALUES,SCATTER_REVERSE); CHKERRQ(ierr);
286  ierr = VecAssemblyBegin(F[ii]); CHKERRQ(ierr);
287  ierr = VecAssemblyEnd(F[ii]); CHKERRQ(ierr);
288  }
289 
290 // ierr = VecView(F,PETSC_VIEWER_STDOUT_WORLD); CHKERRQ(ierr);
291 // ierr = MatView(Aij,PETSC_VIEWER_STDOUT_WORLD); CHKERRQ(ierr);
292 
293 
294  // Volume calculation
295  //==============================================================================================================================
296  VolumeElementForcesAndSourcesCore vol_elem(m_field);
297  Vec volume_vec;
298  int volume_vec_ghost[] = { 0 };
299  ierr = VecCreateGhost(
300  PETSC_COMM_WORLD,(!m_field.get_comm_rank())?1:0,1,1,volume_vec_ghost,&volume_vec
301  ); CHKERRQ(ierr);
302  ierr = VecZeroEntries(volume_vec); CHKERRQ(ierr);
303  vol_elem.getOpPtrVector().push_back(new VolumeCalculation("DISPLACEMENT",volume_vec));
304 
305  ierr = m_field.loop_finite_elements("ELASTIC_MECHANICS","ELASTIC", vol_elem); CHKERRQ(ierr);
306 
307  ierr = VecAssemblyBegin(volume_vec); CHKERRQ(ierr);
308  ierr = VecAssemblyEnd(volume_vec); CHKERRQ(ierr);
309  double rve_volume;
310  ierr = VecSum(volume_vec,&rve_volume); CHKERRQ(ierr);
311 
312  ierr = VecView(volume_vec,PETSC_VIEWER_STDOUT_WORLD); CHKERRQ(ierr);
313  ierr = PetscPrintf(PETSC_COMM_WORLD,"RVE Volume %3.2g\n",rve_volume); CHKERRQ(ierr);
314  ierr = VecDestroy(&volume_vec);
315 
316  //==============================================================================================================================
317  //Post processing
318  //==============================================================================================================================
319  struct MyPostProc: public PostProcVolumeOnRefinedMesh {
320  bool doPreProcess;
321  bool doPostProcess;
322  MyPostProc(MoFEM::Interface &m_field):
323  PostProcVolumeOnRefinedMesh(m_field),
324  doPreProcess(true),
325  doPostProcess(true)
326  {}
327  void setDoPreProcess() { doPreProcess = true; }
328  void unSetDoPreProcess() { doPreProcess = false; }
329  void setDoPostProcess() { doPostProcess = true; }
330  void unSetDoPostProcess() { doPostProcess = false; }
331 
332 
333 
334  PetscErrorCode preProcess() {
335  PetscFunctionBegin;
336  if(doPreProcess) {
337  ierr = PostProcVolumeOnRefinedMesh::preProcess(); CHKERRQ(ierr);
338  }
339  PetscFunctionReturn(0);
340  }
341  PetscErrorCode postProcess() {
342  PetscFunctionBegin;
343  if(doPostProcess) {
344  ierr = PostProcVolumeOnRefinedMesh::postProcess(); CHKERRQ(ierr);
345  }
346  PetscFunctionReturn(0);
347  }
348  };
349 
350  MyPostProc post_proc(m_field);
351 
352  ierr = post_proc.generateReferenceElementMesh(); CHKERRQ(ierr);
353  ierr = post_proc.addFieldValuesPostProc("DISPLACEMENT"); CHKERRQ(ierr);
354  ierr = post_proc.addFieldValuesGradientPostProc("DISPLACEMENT"); CHKERRQ(ierr);
355  ierr = post_proc.addFieldValuesPostProc("MESH_NODE_POSITIONS"); CHKERRQ(ierr);
356  ierr = post_proc.addFieldValuesGradientPostProc("MESH_NODE_POSITIONS"); CHKERRQ(ierr);
357 
358  for(
359  map<int,NonlinearElasticElement::BlockData>::iterator sit = elastic.setOfBlocks.begin();
360  sit != elastic.setOfBlocks.end(); sit++
361  ) {
362  post_proc.getOpPtrVector().push_back(
363  new PostPorcStress(
364  post_proc.postProcMesh,
365  post_proc.mapGaussPts,
366  "DISPLACEMENT",
367  sit->second,
368  post_proc.commonData,
369  false
370  )
371  );
372  }
373  ierr = VecGhostUpdateBegin(D,INSERT_VALUES,SCATTER_FORWARD); CHKERRQ(ierr);
374  ierr = VecGhostUpdateEnd(D,INSERT_VALUES,SCATTER_FORWARD); CHKERRQ(ierr);
375  ierr = m_field.set_global_ghost_vector(
376  "ELASTIC_MECHANICS",ROW,D,INSERT_VALUES,SCATTER_REVERSE
377  ); CHKERRQ(ierr);
378 
379  //==============================================================================================================================
380 
381  //Solver
382  KSP solver;
383  ierr = KSPCreate(PETSC_COMM_WORLD,&solver); CHKERRQ(ierr);
384  ierr = KSPSetOperators(solver,Aij,Aij); CHKERRQ(ierr);
385  ierr = KSPSetFromOptions(solver); CHKERRQ(ierr);
386  ierr = KSPSetUp(solver); CHKERRQ(ierr);
387 
388  MatrixDouble Dmat;
389  Dmat.resize(6,6);
390  Dmat.clear();
391 
392  //calculate honmogenised stress
393  //create a vector for 6 components of homogenized stress
394  Vec stress_homo;
395  int stress_homo_ghost[] = { 0,1,2,3,4,5,6 };
396  ierr = VecCreateGhost(
397  PETSC_COMM_WORLD,(!m_field.get_comm_rank())?6:0,6,6,stress_homo_ghost,&stress_homo
398  ); CHKERRQ(ierr);
399 
400  lagrangian_element.setRVEBCsHomoStressOperators("DISPLACEMENT","LAGRANGE_MUL_DISP","MESH_NODE_POSITIONS",stress_homo);
401 
402  PetscScalar *avec;
403  ierr = VecGetArray(stress_homo,&avec); CHKERRQ(ierr);
404  for(int ii = 0;ii<6;ii++) {
405  ierr = VecZeroEntries(D); CHKERRQ(ierr);
406  ierr = KSPSolve(solver,F[ii],D); CHKERRQ(ierr);
407  ierr = VecGhostUpdateBegin(D,INSERT_VALUES,SCATTER_FORWARD); CHKERRQ(ierr);
408  ierr = VecGhostUpdateEnd(D,INSERT_VALUES,SCATTER_FORWARD); CHKERRQ(ierr);
409  ierr = m_field.set_global_ghost_vector(
410  "ELASTIC_MECHANICS",ROW,D,INSERT_VALUES,SCATTER_REVERSE
411  ); CHKERRQ(ierr);
412 
413  //post-processing the resutls
414  post_proc.setDoPreProcess();
415  post_proc.unSetDoPostProcess();
416  ierr = m_field.loop_finite_elements(
417  "ELASTIC_MECHANICS","ELASTIC",post_proc
418  ); CHKERRQ(ierr);
419  post_proc.unSetDoPreProcess();
420  post_proc.setDoPostProcess();
421  {
422  ostringstream sss;
423  sss << "mode_" << "Disp" << "_" << ii << ".h5m";
424  rval = post_proc.postProcMesh.write_file(
425  sss.str().c_str(),"MOAB","PARALLEL=WRITE_PART"
426  ); CHKERRQ_MOAB(rval);
427  }
428 
429  ierr = VecZeroEntries(stress_homo); CHKERRQ(ierr);
430  ierr = m_field.loop_finite_elements(
431  "ELASTIC_MECHANICS","LAGRANGE_ELEM",lagrangian_element.getLoopFeRVEBCStress()
432  ); CHKERRQ(ierr);
433  ierr = PetscOptionsGetReal(
434  PETSC_NULL,"-my_rve_volume",&rve_volume,PETSC_NULL
435  ); CHKERRQ(ierr);
436  ierr = VecAssemblyBegin(stress_homo); CHKERRQ(ierr);
437  ierr = VecAssemblyEnd(stress_homo); CHKERRQ(ierr);
438  ierr = VecScale(stress_homo,1.0/rve_volume); CHKERRQ(ierr);
439  ierr = VecGhostUpdateBegin(stress_homo,INSERT_VALUES,SCATTER_FORWARD); CHKERRQ(ierr);
440  ierr = VecGhostUpdateEnd(stress_homo,INSERT_VALUES,SCATTER_FORWARD); CHKERRQ(ierr);
441  for(int jj=0; jj<6; jj++) {
442  Dmat(jj,ii) = avec[jj];
443  }
444  }
445 
446  ierr = VecRestoreArray(stress_homo,&avec); CHKERRQ(ierr);
447  PetscPrintf(
448  PETSC_COMM_WORLD,"\nHomogenised Stiffens Matrix = \n\n"
449  );
450 
451  for(int ii=0; ii<6; ii++) {
452  PetscPrintf(
453  PETSC_COMM_WORLD,
454  "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",
455  ii,Dmat(ii,0),Dmat(ii,1),Dmat(ii,2),Dmat(ii,3),Dmat(ii,4),Dmat(ii,5)
456  );
457  }
458  //==============================================================================================================================
459 
460 
461 
462  //Open mesh_file_name.txt for writing
463  ofstream myfile;
464  myfile.open ((string(mesh_file_name)+".txt").c_str());
465 
466  //Output displacements
467  myfile << "<<<< Homonenised stress >>>>>" << endl;
468 
469  if(pcomm->rank()==0){
470  for(int ii=0; ii<6; ii++){
471  myfile << boost::format("%.3lf") % roundn(Dmat(ii,0)) << endl;
472  avec++;
473  }
474  }
475  //Close mesh_file_name.txt
476  myfile.close();
477 
478 
479 
480 
481 
482  //detroy matrices
483  for(int ii = 0;ii<6;ii++) {
484  ierr = VecDestroy(&F[ii]); CHKERRQ(ierr);
485  }
486  ierr = VecDestroy(&D); CHKERRQ(ierr);
487  ierr = MatDestroy(&Aij); CHKERRQ(ierr);
488  ierr = KSPDestroy(&solver); CHKERRQ(ierr);
489  ierr = VecDestroy(&stress_homo); CHKERRQ(ierr);
490 
491  ierr = PetscTime(&v2);CHKERRQ(ierr);
492  ierr = PetscGetCPUTime(&t2);CHKERRQ(ierr);
493 
494  PetscSynchronizedPrintf(PETSC_COMM_WORLD,"Total Rank %d Time = %f CPU Time = %f\n",pcomm->rank(),v2-v1,t2-t1);
495 
496  }
497  CATCH_ERRORS;
498 
499 
500  MoFEM::Core::Finalize();
501 
502 }
MethodForForceScaling.hpp
SIDESET
@ SIDESET
Definition: definitions.h:147
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:215
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:454
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
RND_EPS
#define RND_EPS
Definition: homogenisation_disp_atom_test.cpp:48
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
BCs_RVELagrange_Disp::setRVEBCsOperators
PetscErrorCode setRVEBCsOperators(string field_name, string lagrang_field_name, string mesh_nodals_positions, Mat aij, vector< Vec > &fvec)
Definition: BCs_RVELagrange_Disp.hpp:580
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:123
MoFEM::Interface
DeprecatedCoreInterface Interface
Definition: Interface.hpp:1975
BCs_RVELagrange_Disp
Definition: BCs_RVELagrange_Disp.hpp:18
MoFEM::Exceptions::rval
static MoFEMErrorCodeGeneric< moab::ErrorCode > rval
Definition: Exceptions.hpp:74
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
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.
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.
TimeForceScale.hpp
roundn
double roundn(double n)
Definition: homogenisation_disp_atom_test.cpp:52
COL
@ COL
Definition: definitions.h:123
NonLinearElasticElement.hpp
Operators and data structures for non-linear elastic analysis.
main
int main(int argc, char *argv[])
Definition: homogenisation_disp_atom_test.cpp:84
Projection10NodeCoordsOnField.hpp
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.
help
static char help[]
Definition: homogenisation_disp_atom_test.cpp:47
MyPostProc::generateReferenceElementMesh
MoFEMErrorCode generateReferenceElementMesh()
Definition: plot_base.cpp:397
MoFEM::VolumeElementForcesAndSourcesCore
Volume finite element base.
Definition: VolumeElementForcesAndSourcesCore.hpp:26
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
mesh_file_name
char mesh_file_name[255]
Definition: mesh_smoothing.cpp:23
convert.n
n
Definition: convert.py:82
Range
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
CATCH_ERRORS
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:372
PostProcVolumeOnRefinedMesh
Post processing.
Definition: PostProcOnRefMesh.hpp:955
PostPorcStress
DEPRECATED typedef PostProcStress PostPorcStress
Definition: PostProcStresses.hpp:193
MoFEM::Core
CoreTmp< 0 > Core
Definition: Core.hpp:1094
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
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
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
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
Generated by Doxygen 1.8.17 and hosted at University of Glasgow