84 {
85
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
96 PetscBool flg = PETSC_TRUE;
99 if(flg != PETSC_TRUE) {
100 SETERRQ(PETSC_COMM_SELF,1,"*** ERROR -my_file (MESH FILE NEEDED)");
101 }
104 if(flg != PETSC_TRUE) {
106 }
107
108
109 const char *option;
110 option = "";
112 ParallelComm* pcomm = ParallelComm::get_pcomm(&moab,
MYPCOMM_INDEX);
113 if(pcomm == NULL) pcomm = new ParallelComm(&moab,PETSC_COMM_WORLD);
114
115
116 PetscLogDouble t1,t2;
117 PetscLogDouble v1,v2;
118 ierr = PetscTime(&v1); CHKERRQ(
ierr);
119 ierr = PetscGetCPUTime(&t1); CHKERRQ(
ierr);
120
121
124
125
127
128
130 bit_level0.set(0);
134 ierr = m_field.get_entities_by_ref_level(bit_level0,
BitRefLevel().set(),meshset_level0); CHKERRQ(
ierr);
135
136
137
138
139
140 int field_rank=3;
144
145
146
149
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
157 ierr = m_field.seed_ref_level_MESHSET(BoundFacesMeshset,
BitRefLevel().set()); CHKERRQ(
ierr);
159
160
161
162
167
172
176
177
178
179
180 boost::shared_ptr<Hooke<adouble> > hooke_adouble(
new Hooke<adouble>);
181 boost::shared_ptr<Hooke<double> > hooke_double(
new Hooke<double>);
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
189 lagrangian_element.addLagrangiangElement("LAGRANGE_ELEM","DISPLACEMENT","LAGRANGE_MUL_DISP","MESH_NODE_POSITIONS");
190
191
193
194
197
198
199
201
202
203
204
205
206
207
208
211 ierr = m_field.
loop_dofs(
"MESH_NODE_POSITIONS",ent_method_material); CHKERRQ(
ierr);
212
213
215
216
218
219
220 ierr = m_field.build_problems(); CHKERRQ(
ierr);
221
222
223
224
225
226
227 ierr = m_field.partition_problem(
"ELASTIC_MECHANICS"); CHKERRQ(
ierr);
228 ierr = m_field.partition_finite_elements(
"ELASTIC_MECHANICS"); CHKERRQ(
ierr);
229
230 ierr = m_field.partition_ghost_dofs(
"ELASTIC_MECHANICS"); CHKERRQ(
ierr);
231
232
233 ierr = m_field.print_cubit_displacement_set(); CHKERRQ(
ierr);
234 ierr = m_field.print_cubit_force_set(); CHKERRQ(
ierr);
235
236 ierr = m_field.print_cubit_materials_set(); CHKERRQ(
ierr);
237
238
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
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
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
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
266 elastic.getLoopFeLhs().snes_B = Aij;
268
269 lagrangian_element.setRVEBCsOperators(
"DISPLACEMENT",
"LAGRANGE_MUL_DISP",
"MESH_NODE_POSITIONS",Aij,
F);
271 cout<<"after lagrangian_element.getLoopFeRVEBCLhs "<<endl;
273 cout<<"after lagrangian_element.getLoopFeRVEBCRhs "<<endl;
274
275
276
277
278
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
291
292
293
294
295
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
302 ierr = VecZeroEntries(volume_vec); CHKERRQ(
ierr);
303 vol_elem.getOpPtrVector().push_back(
new VolumeCalculation(
"DISPLACEMENT",volume_vec));
304
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
318
320 bool doPreProcess;
321 bool doPostProcess;
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
335 PetscFunctionBegin;
336 if(doPreProcess) {
338 }
339 PetscFunctionReturn(0);
340 }
342 PetscFunctionBegin;
343 if(doPostProcess) {
345 }
346 PetscFunctionReturn(0);
347 }
348 };
349
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(
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
378
379
380
381
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
389 Dmat.resize(6,6);
390 Dmat.clear();
391
392
393
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
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++) {
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
412
413
414 post_proc.setDoPreProcess();
415 post_proc.unSetDoPostProcess();
417 "ELASTIC_MECHANICS","ELASTIC",post_proc
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"
427 }
428
429 ierr = VecZeroEntries(stress_homo); CHKERRQ(
ierr);
431 "ELASTIC_MECHANICS","LAGRANGE_ELEM",lagrangian_element.getLoopFeRVEBCStress()
434 PETSC_NULL,"-my_rve_volume",&rve_volume,PETSC_NULL
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
463 ofstream myfile;
465
466
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
476 myfile.close();
477
478
479
480
481
482
483 for(int ii = 0;ii<6;ii++) {
484 ierr = VecDestroy(&
F[ii]); CHKERRQ(
ierr);
485 }
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 }
498
499
501
502}
DEPRECATED typedef PostProcStress PostPorcStress
#define CATCH_ERRORS
Catch errors.
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
#define MYPCOMM_INDEX
default communicator number PCOMM
#define CHKERRQ_MOAB(a)
check error code of MoAB function
constexpr int order
Order displacement.
virtual MoFEMErrorCode build_finite_elements(int verb=DEFAULT_VERBOSITY)=0
Build finite elements.
virtual MoFEMErrorCode build_fields(int verb=DEFAULT_VERBOSITY)=0
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.
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.
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.
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.
virtual MoFEMErrorCode add_problem(const std::string &name, enum MoFEMTypes bh=MF_EXCL, int verb=DEFAULT_VERBOSITY)=0
Add problem.
virtual MoFEMErrorCode modify_problem_ref_level_add_bit(const std::string &name_problem, const BitRefLevel &bit)=0
add ref level to problem
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
const FTensor::Tensor2< T, Dim, Dim > Vec
static MoFEMErrorCodeGeneric< PetscErrorCode > ierr
static MoFEMErrorCodeGeneric< moab::ErrorCode > rval
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
PetscErrorCode PetscOptionsGetInt(PetscOptions *, const char pre[], const char name[], PetscInt *ivalue, PetscBool *set)
PetscErrorCode PetscOptionsGetReal(PetscOptions *, const char pre[], const char name[], PetscReal *dval, PetscBool *set)
PetscErrorCode PetscOptionsGetString(PetscOptions *, const char pre[], const char name[], char str[], size_t size, PetscBool *set)
virtual MoFEMErrorCode build_adjacencies(const Range &ents, int verb=DEFAULT_VERBOSITY)=0
build adjacencies
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.
virtual int get_comm_rank() const =0
static MoFEMErrorCode Initialize(int *argc, char ***args, const char file[], const char help[])
Initializes the MoFEM database PETSc, MOAB and MPI.
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Deprecated interface functions.
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...
Projection of edge entities with one mid-node on hierarchical basis.
Volume finite element base.
MoFEMErrorCode postProcess()
MoFEMErrorCode preProcess()
structure grouping operators and data used for calculation of nonlinear elastic element
MoFEMErrorCode preProcess()
MoFEMErrorCode postProcess()