v0.13.1
Functions | Variables
rve_mech_plas_interface_disp.cpp File Reference
#include <MoFEM.hpp>
#include <boost/program_options.hpp>
#include <boost/numeric/ublas/vector_proxy.hpp>
#include <boost/numeric/ublas/matrix.hpp>
#include <boost/numeric/ublas/matrix_proxy.hpp>
#include <boost/numeric/ublas/vector.hpp>
#include <boost/numeric/ublas/symmetric.hpp>
#include <MethodForForceScaling.hpp>
#include <TimeForceScale.hpp>
#include <BCs_RVELagrange_Disp.hpp>
#include <VolumeCalculation.hpp>
#include <CohesiveInterfaceElement.hpp>
#include <adolc/adolc.h>
#include <NonLinearElasticElement.hpp>
#include <SmallTransverselyIsotropic.hpp>
#include <SmallStrainPlasticity.hpp>
#include <SmallStrainPlasticityMaterialModels.hpp>
#include <Projection10NodeCoordsOnField.hpp>
#include <PostProcOnRefMesh.hpp>
#include <string>

Go to the source code of this file.

Functions

int main (int argc, char *argv[])
 

Variables

static char help []
 

Function Documentation

◆ main()

int main ( int  argc,
char *  argv[] 
)

Definition at line 70 of file rve_mech_plas_interface_disp.cpp.

70 {
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;
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);
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
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;
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;
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();
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
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(
584 "DISPLACEMENT",small_strain_plasticity.commonData
585 )
586 );
587 small_strain_plasticity.feRhs.getOpPtrVector().push_back(
589 m_field,"DISPLACEMENT",small_strain_plasticity.commonData,cp
590 )
591 );
592 small_strain_plasticity.feRhs.getOpPtrVector().push_back(
594 "DISPLACEMENT",small_strain_plasticity.commonData
595 )
596 );
597 small_strain_plasticity.feLhs.getOpPtrVector().push_back(
599 "DISPLACEMENT",small_strain_plasticity.commonData
600 )
601 );
602 small_strain_plasticity.feLhs.getOpPtrVector().push_back(
604 m_field,"DISPLACEMENT",small_strain_plasticity.commonData,cp
605 )
606 );
607 small_strain_plasticity.feLhs.getOpPtrVector().push_back(
609 "DISPLACEMENT",small_strain_plasticity.commonData
610 )
611 );
612 small_strain_plasticity.feUpdate.getOpPtrVector().push_back(
614 "DISPLACEMENT",small_strain_plasticity.commonData
615 )
616 );
617 small_strain_plasticity.feUpdate.getOpPtrVector().push_back(
619 m_field,"DISPLACEMENT",small_strain_plasticity.commonData,cp
620 )
621 );
622 small_strain_plasticity.feUpdate.getOpPtrVector().push_back(
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
765 dm,D,INSERT_VALUES,SCATTER_REVERSE
766 ); CHKERRQ(ierr);
767
769 dm,"PLASTIC",&small_strain_plasticity.feUpdate
770 ); CHKERRQ(ierr);
771
772
773
775 dm,"INTERFACE",&cohesive_elements.feHistory
776 ); CHKERRQ(ierr);
777
778
779 //Save data on mesh
780 {
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);
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 }
871
873}
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:385
@ MF_ZERO
Definition: definitions.h:111
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:73
@ H1
continuous field
Definition: definitions.h:98
#define MYPCOMM_INDEX
default communicator number PCOMM
Definition: definitions.h:228
@ SIDESET
Definition: definitions.h:160
@ UNKNOWNNAME
Definition: definitions.h:171
@ BLOCKSET
Definition: definitions.h:161
#define CHKERRQ_MOAB(a)
check error code of MoAB function
Definition: definitions.h:467
#define MU(E, NU)
Definition: fem_tools.h:33
#define LAMBDA(E, NU)
Definition: fem_tools.h:32
PetscErrorCode DMMoFEMSetIsPartitioned(DM dm, PetscBool is_partitioned)
Definition: DMMMoFEM.cpp:1082
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: DMMMoFEM.cpp:677
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: DMMMoFEM.cpp:130
PetscErrorCode DMMoFEMAddElement(DM dm, const char fe_name[])
add element to dm
Definition: DMMMoFEM.cpp:462
PetscErrorCode DMoFEMMeshToLocalVector(DM dm, Vec l, InsertMode mode, ScatterMode scatter_mode)
set local (or ghosted) vector values on mesh for partition only
Definition: DMMMoFEM.cpp:482
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: DMMMoFEM.cpp:718
PetscErrorCode DMCreateMatrix_MoFEM(DM dm, Mat *M)
Definition: DMMMoFEM.cpp:1156
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition: DMMMoFEM.cpp:59
PetscErrorCode DMoFEMLoopFiniteElements(DM dm, const char fe_name[], MoFEM::FEMethod *method, CacheTupleWeakPtr cache_ptr=CacheTupleSharedPtr())
Executes FEMethod for finite elements in DM.
Definition: DMMMoFEM.cpp:545
PetscErrorCode DMCreateGlobalVector_MoFEM(DM dm, Vec *g)
DMShellSetCreateGlobalVector.
Definition: DMMMoFEM.cpp:1126
PetscErrorCode DMMoFEMGetSnesCtx(DM dm, MoFEM::SnesCtx **snes_ctx)
get MoFEM::SnesCtx data structure
Definition: DMMMoFEM.cpp:1053
PetscErrorCode DMoFEMMeshToGlobalVector(DM dm, Vec g, InsertMode mode, ScatterMode scatter_mode)
set ghosted vector values on all existing mesh entities
Definition: DMMMoFEM.cpp:494
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
virtual MoFEMErrorCode add_finite_element(const std::string &fe_name, enum MoFEMTypes bh=MF_EXCL, int verb=DEFAULT_VERBOSITY)=0
add finite element
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
virtual MoFEMErrorCode build_finite_elements(int verb=DEFAULT_VERBOSITY)=0
Build finite elements.
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
virtual MoFEMErrorCode modify_finite_element_add_field_data(const std::string &fe_name, const std::string &name_filed)=0
set finite element field data
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
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 bool check_field(const std::string &name) const =0
check if field is in database
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.
#define _IT_CUBITMESHSETS_BY_BCDATA_TYPE_FOR_LOOP_(MESHSET_MANAGER, CUBITBCTYPE, IT)
Iterator that loops over a specific Cubit MeshSet in a moFEM field.
#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.
char mesh_file_name[255]
const FTensor::Tensor2< T, Dim, Dim > Vec
static MoFEMErrorCodeGeneric< PetscErrorCode > ierr
Definition: Exceptions.hpp:87
static MoFEMErrorCodeGeneric< moab::ErrorCode > rval
Definition: Exceptions.hpp:85
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
Definition: Types.hpp:51
UBlasVector< double > VectorDouble
Definition: Types.hpp:79
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:148
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 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:39
PetscErrorCode PetscOptionsGetBool(PetscOptions *, const char pre[], const char name[], PetscBool *bval, PetscBool *set)
PetscErrorCode PetscOptionsGetRealArray(PetscOptions *, const char pre[], const char name[], PetscReal dval[], PetscInt *nmax, PetscBool *set)
CoreTmp< 0 > Core
Definition: Core.hpp:1096
PetscErrorCode PetscOptionsGetString(PetscOptions *, const char pre[], const char name[], char str[], size_t size, PetscBool *set)
DeprecatedCoreInterface Interface
Definition: Interface.hpp:1965
const double D
diffusivity
double young_modulus
Definition: plastic.cpp:105
double scale
Definition: plastic.cpp:103
constexpr double t
plate stiffness
Definition: plate.cpp:76
static char help[]
virtual moab::Interface & get_moab()=0
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
Core (interface) class.
Definition: Core.hpp:92
static MoFEMErrorCode Initialize(int *argc, char ***args, const char file[], const char help[])
Initializes the MoFEM database PETSc, MOAB and MPI.
Definition: Core.cpp:85
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition: Core.cpp:125
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...
Transverse Isotropic material data structure.
Linear interface data structure.
Projection of edge entities with one mid-node on hierarchical basis.
Interface for nonlinear (SNES) solver.
Definition: SnesCtx.hpp:27
structure grouping operators and data used for calculation of nonlinear elastic element
Post processing.
J2 plasticity (Kinematic Isotropic (Linear) Hardening)
Small strain finite element implementation.
Force scale operator for reading two columns.
Calculate volume.

Variable Documentation

◆ help

char help[]
static
Initial value:
=
"...\n"
"\n"

Definition at line 66 of file rve_mech_plas_interface_disp.cpp.