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

Variable Documentation

◆ help

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

Definition at line 66 of file rve_mech_plas_interface_disp.cpp.

Generated by Doxygen 1.12.0 and hosted at University of Glasgow