v0.15.0
Loading...
Searching...
No Matches
ElasticMaterials.hpp
Go to the documentation of this file.
1/** \file ElasticMaterials.hpp
2 * \ingroup nonlinear_elastic_elem
3 * \brief Elastic materials
4 */
5
6
7
8#ifndef __ELASTICMATERIALS_HPP__
9#define __ELASTICMATERIALS_HPP__
10
11#include <boost/program_options.hpp>
12namespace po = boost::program_options;
13
14#include <Hooke.hpp>
15#include <NeoHookean.hpp>
16
17#define MAT_KIRCHHOFF "KIRCHHOFF"
18#define MAT_HOOKE "HOOKE"
19#define MAT_NEOHOOKEAN "NEOHOOKEAN"
20
21/** \brief Manage setting parameters and constitutive equations for
22 * nonlinear/linear elastic materials \ingroup nonlinear_elastic_elem
23 */
25
27 std::string defMaterial; ///< default material, if block is set to elastic,
28 ///< this material is used as default
29 std::string configFile; //< default name of config file
30
31 bool iNitialized; ///< true if class is initialized
32
34 std::string def_material = MAT_KIRCHHOFF)
35 : mField(m_field), defMaterial(def_material),
36 configFile("elastic_material.in"), iNitialized(false) {}
37
38 virtual ~ElasticMaterials() {}
39
40 std::map<std::string,
41 boost::shared_ptr<NonlinearElasticElement::
42 FunctionsToCalculatePiolaKirchhoffI<adouble>>>
43 aDoubleMaterialModel; ///< Hash map of materials for evaluation with
44 ///< adouble, i.e. ADOL-C
45
46 std::map<
47 std::string,
48 boost::shared_ptr<
50 doubleMaterialModel; ///< Hash map of materials for evaluation with double
51
52 /**
53 * Structure for material parameters for block
54 */
56 string mAterial;
57 int oRder;
58 double yOung;
59 double pOisson;
60 double dEnsity;
61 double dashG;
63 double aX, aY, aZ;
65 : mAterial(MAT_KIRCHHOFF), oRder(-1), yOung(-1), pOisson(-2),
66 dEnsity(-1), dashG(-1), dashPoisson(-1), aX(0), aY(0), aZ(0) {}
67 };
68
69 std::map<int, BlockOptionData> blockData; ///< Material parameters on blocks
70
71 PetscBool isConfigFileSet; ///< True if config file is set from command line
72 po::variables_map vM;
73
74 /**
75 * Initialize model parameters
76 * @return [description]
77 */
78 virtual MoFEMErrorCode iNit() {
80 // add new material below
81 string mat_name;
83 boost::make_shared<NonlinearElasticElement::
84 FunctionsToCalculatePiolaKirchhoffI<adouble>>();
85 doubleMaterialModel[MAT_KIRCHHOFF] = boost::make_shared<
87 aDoubleMaterialModel[MAT_HOOKE] = boost::make_shared<Hooke<adouble>>();
88 doubleMaterialModel[MAT_HOOKE] = boost::make_shared<Hooke<double>>();
90 boost::make_shared<NeoHookean<adouble>>();
92 boost::make_shared<NeoHookean<double>>();
93 std::ostringstream avilable_materials;
94 avilable_materials << "set elastic material < ";
95 std::map<std::string,
96 boost::shared_ptr<
98 double>>>::iterator mit;
99 mit = doubleMaterialModel.begin();
100 for (; mit != doubleMaterialModel.end(); mit++) {
101 avilable_materials << mit->first << " ";
102 }
103 avilable_materials << ">";
104
105 PetscOptionsBegin(mField.get_comm(), "", "Elastic Materials Configuration",
106 "none");
107 char default_material[255];
108 PetscBool def_mat_set;
109 CHKERR PetscOptionsString(
110 "-default_material", avilable_materials.str().c_str(), "",
111 defMaterial.c_str(), default_material, 255, &def_mat_set);
112 if (def_mat_set) {
113 defMaterial = default_material;
115 aDoubleMaterialModel.end()) {
116 SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
117 "material <%s> not implemented", default_material);
118 }
119 }
120 char config_file[255];
121 CHKERR PetscOptionsString("-elastic_material_configuration",
122 "elastic materials configure file name", "",
123 configFile.c_str(), config_file, 255,
125 if (isConfigFileSet) {
126 configFile = config_file;
127 }
128 PetscOptionsEnd();
130 }
131
132 /** \brief read Elastic materials declaration for blocks and meshsets
133
134 File parameters:
135 \code
136 [block_1]
137 displacemet_order = 1/2 .. N
138 material = KIRCHHOFF/HOOKE/NEOHOOKEAN
139 young_modulus = 1
140 poisson_ratio = 0.25
141 density = 1
142 a_x = 0
143 a_y = 0
144 a_z = 10
145 \endcode
146
147 To read material configuration file you need to use option:
148 \code
149 -elastic_material_configuration name_of_config_file
150 \endcode
151
152 */
153 virtual MoFEMErrorCode readConfigFile() {
155
156 ifstream file(configFile.c_str());
157 if (isConfigFileSet) {
158 if (!file.good()) {
159 SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_FOUND, "file < %s > not found",
160 configFile.c_str());
161 }
162 } else {
163 if (!file.good()) {
165 }
166 }
167
168 po::options_description config_file_options;
170
171 std::ostringstream str_order;
172 str_order << "block_" << it->getMeshsetId() << ".displacemet_order";
173 config_file_options.add_options()(
174 str_order.str().c_str(),
175 po::value<int>(&blockData[it->getMeshsetId()].oRder)
176 ->default_value(-1));
177
178 std::ostringstream str_material;
179 str_material << "block_" << it->getMeshsetId() << ".material";
180 config_file_options.add_options()(
181 str_material.str().c_str(),
182 po::value<std::string>(&blockData[it->getMeshsetId()].mAterial)
183 ->default_value(defMaterial));
184
185 std::ostringstream str_ym;
186 str_ym << "block_" << it->getMeshsetId() << ".young_modulus";
187 config_file_options.add_options()(
188 str_ym.str().c_str(),
189 po::value<double>(&blockData[it->getMeshsetId()].yOung)
190 ->default_value(-1));
191
192 std::ostringstream str_pr;
193 str_pr << "block_" << it->getMeshsetId() << ".poisson_ratio";
194 config_file_options.add_options()(
195 str_pr.str().c_str(),
196 po::value<double>(&blockData[it->getMeshsetId()].pOisson)
197 ->default_value(-2));
198
199 std::ostringstream str_density;
200 str_density << "block_" << it->getMeshsetId() << ".density";
201 config_file_options.add_options()(
202 str_density.str().c_str(),
203 po::value<double>(&blockData[it->getMeshsetId()].dEnsity)
204 ->default_value(-1));
205
206 std::ostringstream str_dashG;
207 str_dashG << "block_" << it->getMeshsetId() << ".dashG";
208 config_file_options.add_options()(
209 str_dashG.str().c_str(),
210 po::value<double>(&blockData[it->getMeshsetId()].dashG)
211 ->default_value(-1));
212
213 std::ostringstream str_dashPoisson;
214 str_dashPoisson << "block_" << it->getMeshsetId() << ".dashPoisson";
215 config_file_options.add_options()(
216 str_dashPoisson.str().c_str(),
217 po::value<double>(&blockData[it->getMeshsetId()].dashPoisson)
218 ->default_value(-2));
219
220 std::ostringstream str_ax;
221 str_ax << "block_" << it->getMeshsetId() << ".a_x";
222 config_file_options.add_options()(
223 str_ax.str().c_str(),
224 po::value<double>(&blockData[it->getMeshsetId()].aX)
225 ->default_value(0));
226
227 std::ostringstream str_ay;
228 str_ay << "block_" << it->getMeshsetId() << ".a_y";
229 config_file_options.add_options()(
230 str_ay.str().c_str(),
231 po::value<double>(&blockData[it->getMeshsetId()].aY)
232 ->default_value(0));
233
234 std::ostringstream str_az;
235 str_az << "block_" << it->getMeshsetId() << ".a_z";
236 config_file_options.add_options()(
237 str_az.str().c_str(),
238 po::value<double>(&blockData[it->getMeshsetId()].aZ)
239 ->default_value(0));
240 }
241 po::parsed_options parsed =
242 parse_config_file(file, config_file_options, true);
243 store(parsed, vM);
244 po::notify(vM);
245 std::vector<std::string> additional_parameters;
246 additional_parameters =
247 collect_unrecognized(parsed.options, po::include_positional);
248 for (std::vector<std::string>::iterator vit = additional_parameters.begin();
249 vit != additional_parameters.end(); vit++) {
250 CHKERR PetscPrintf(PETSC_COMM_WORLD,
251 "** WARNING Unrecognized option %s\n", vit->c_str());
252 }
253
255 }
256
257 MoFEMErrorCode setBlocksOrder() {
259
260 // set app. order
261 PetscBool flg = PETSC_TRUE;
262 PetscInt disp_order;
263 CHKERR PetscOptionsGetInt(PETSC_NULLPTR, PETSC_NULLPTR, "-order", &disp_order,
264 &flg);
265 if (flg != PETSC_TRUE) {
266 disp_order = 1;
267 }
269 if (blockData[it->getMeshsetId()].oRder == -1)
270 continue;
271 if (blockData[it->getMeshsetId()].oRder == disp_order)
272 continue;
273 PetscPrintf(mField.get_comm(), "Set block %d oRder to %d\n",
274 it->getMeshsetId(), blockData[it->getMeshsetId()].oRder);
275 Range block_ents;
276 rval = mField.get_moab().get_entities_by_handle(it->meshset, block_ents,
277 true);
278 CHKERRG(rval);
279 Range ents_to_set_order;
280 CHKERR mField.get_moab().get_adjacencies(
281 block_ents, 3, false, ents_to_set_order, moab::Interface::UNION);
282 ents_to_set_order = ents_to_set_order.subset_by_type(MBTET);
283 CHKERR mField.get_moab().get_adjacencies(
284 block_ents, 2, false, ents_to_set_order, moab::Interface::UNION);
285 CHKERR mField.get_moab().get_adjacencies(
286 block_ents, 1, false, ents_to_set_order, moab::Interface::UNION);
287 if (mField.check_field("DISPLACEMENT")) {
288 CHKERR mField.set_field_order(ents_to_set_order, "DISPLACEMENT",
289 blockData[it->getMeshsetId()].oRder);
290 }
291 if (mField.check_field("SPATIAL_POSITION")) {
292 CHKERR mField.set_field_order(ents_to_set_order, "SPATIAL_POSITION",
293 blockData[it->getMeshsetId()].oRder);
294 }
295 if (mField.check_field("DOT_SPATIAL_POSITION")) {
296 CHKERR mField.set_field_order(ents_to_set_order, "DOT_SPATIAL_POSITION",
297 blockData[it->getMeshsetId()].oRder);
298 }
299 if (mField.check_field("SPATIAL_VELOCITY")) {
300 CHKERR mField.set_field_order(ents_to_set_order, "SPATIAL_VELOCITY",
301 blockData[it->getMeshsetId()].oRder);
302 }
303 }
305 }
306
307#ifdef __NONLINEAR_ELASTIC_HPP
308
309 virtual MoFEMErrorCode
310 setBlocks(std::map<int, NonlinearElasticElement::BlockData> &set_of_blocks) {
312
313 if (!iNitialized) {
314 CHKERR iNit();
317 iNitialized = true;
318 }
321 int id = it->getMeshsetId();
322 Mat_Elastic mydata;
323 CHKERR it->getAttributeDataStructure(mydata);
324 EntityHandle meshset = it->getMeshset();
325 CHKERR mField.get_moab().get_entities_by_type(
326 meshset, MBTET, set_of_blocks[id].tEts, true);
327 set_of_blocks[id].iD = id;
328 set_of_blocks[id].E = mydata.data.Young;
329 set_of_blocks[id].PoissonRatio = mydata.data.Poisson;
330 if (blockData[id].yOung >= 0)
331 set_of_blocks[id].E = blockData[id].yOung;
332 if (blockData[id].pOisson >= -1)
333 set_of_blocks[id].PoissonRatio = blockData[id].pOisson;
334
335 blockData[id].mAterial = defMaterial;
336 set_of_blocks[id].materialDoublePtr = doubleMaterialModel.at(defMaterial);
337 set_of_blocks[id].materialAdoublePtr =
339
340 MOFEM_LOG_C("WORLD", Sev::verbose,
341 "Block Id %d Young Modulus %3.2g Poisson Ratio %3.2f "
342 "Material model %s Nb. of elements %d\n",
343 id, set_of_blocks[id].E, set_of_blocks[id].PoissonRatio,
344 blockData[id].mAterial.c_str(),
345 set_of_blocks[id].tEts.size());
346
347 if (blockData[id].mAterial.compare(MAT_KIRCHHOFF) == 0) {
348 set_of_blocks[id].materialDoublePtr =
350 set_of_blocks[id].materialAdoublePtr =
352 } else if (blockData[id].mAterial.compare(MAT_HOOKE) == 0) {
353 set_of_blocks[id].materialDoublePtr = doubleMaterialModel.at(MAT_HOOKE);
354 set_of_blocks[id].materialAdoublePtr =
356 } else if (blockData[id].mAterial.compare(MAT_NEOHOOKEAN) == 0) {
357 set_of_blocks[id].materialDoublePtr =
359 set_of_blocks[id].materialAdoublePtr =
361 } else {
362 SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
363 "field with that space is not implemented");
364 }
365 }
367 }
368
369#endif //__NONLINEAR_ELASTIC_HPP
370
371#ifdef __CONVECTIVE_MASS_ELEMENT_HPP
372
374 setBlocks(std::map<int, ConvectiveMassElement::BlockData> &set_of_blocks) {
376
377 if (!iNitialized) {
378 CHKERR iNit();
381 iNitialized = true;
382 }
384 mField, BLOCKSET | BODYFORCESSET, it)) {
385 int id = it->getMeshsetId();
386 EntityHandle meshset = it->getMeshset();
387 CHKERR mField.get_moab().get_entities_by_type(
388 meshset, MBTET, set_of_blocks[id].tEts, true);
389 Block_BodyForces mydata;
390 CHKERR it->getAttributeDataStructure(mydata);
391 set_of_blocks[id].rho0 = mydata.data.density;
392 set_of_blocks[id].a0.resize(3);
393 set_of_blocks[id].a0[0] = mydata.data.acceleration_x;
394 set_of_blocks[id].a0[1] = mydata.data.acceleration_y;
395 set_of_blocks[id].a0[2] = mydata.data.acceleration_z;
396 if (blockData[id].dEnsity >= 0) {
397 set_of_blocks[id].rho0 = blockData[id].dEnsity;
398 std::ostringstream str_ax;
399 str_ax << "block_" << it->getMeshsetId() << ".a_x";
400 std::ostringstream str_ay;
401 str_ay << "block_" << it->getMeshsetId() << ".a_y";
402 std::ostringstream str_az;
403 str_az << "block_" << it->getMeshsetId() << ".a_z";
404 if (vM.count(str_ax.str().c_str())) {
405 set_of_blocks[id].a0[0] = blockData[id].aX;
406 }
407 if (vM.count(str_ay.str().c_str())) {
408 set_of_blocks[id].a0[1] = blockData[id].aY;
409 }
410 if (vM.count(str_az.str().c_str())) {
411 set_of_blocks[id].a0[2] = blockData[id].aZ;
412 }
413 }
414 PetscPrintf(mField.get_comm(),
415 "Block Id %d Density %3.2g a_x = %3.2g a_y = %3.2g a_z = "
416 "%3.2g Nb. of elements %ld\n",
417 id, set_of_blocks[id].rho0, set_of_blocks[id].a0[0],
418 set_of_blocks[id].a0[1], set_of_blocks[id].a0[2],
419 set_of_blocks[id].tEts.size());
420 }
421
423 }
424
425#endif //__CONVECTIVE_MASS_ELEMENT_HPP
426
427#ifdef __KELVIN_VOIGT_DAMPER_HPP__
428
430 std::map<int, KelvinVoigtDamper::BlockMaterialData> &set_of_blocks) {
432
433 if (!iNitialized) {
434 CHKERR iNit();
437 iNitialized = true;
438 }
439
441 bool set = false;
442 int id = it->getMeshsetId();
443 EntityHandle meshset = it->getMeshset();
444 if (it->getName().compare(0, 6, "DAMPER") == 0) {
445 set = true;
446 std::vector<double> data;
447 CHKERR it->getAttributes(data);
448 if (data.size() < 2) {
449 SETERRQ(PETSC_COMM_SELF, 1, "Data inconsistency");
450 }
451 rval = mField.get_moab().get_entities_by_type(
452 it->meshset, MBTET, set_of_blocks[it->getMeshsetId()].tEts, true);
453 CHKERRG(rval);
454 set_of_blocks[it->getMeshsetId()].gBeta = data[0];
455 set_of_blocks[it->getMeshsetId()].vBeta = data[1];
456 }
457 if (blockData[id].dashG > 0) {
458 set = true;
459 Range tEts;
460 rval =
461 mField.get_moab().get_entities_by_type(meshset, MBTET, tEts, true);
462 CHKERRG(rval);
463 if (tEts.empty())
464 continue;
465 set_of_blocks[it->getMeshsetId()].tEts = tEts;
466 set_of_blocks[it->getMeshsetId()].gBeta = blockData[id].dashG;
467 set_of_blocks[it->getMeshsetId()].vBeta = blockData[id].dashPoisson;
468 }
469 if (set) {
470 PetscPrintf(
472 "Block Id %d Damper Shear Modulus = %3.2g Poisson ratio = %3.2g\n",
473 id, set_of_blocks[id].gBeta, set_of_blocks[id].vBeta);
474 }
475 }
476
478 }
479
480#endif //__KELVIN_VOIGT_DAMPER_HPP__
481};
482
483#endif //__ELASTICMATERIALS_HPP__
#define MAT_KIRCHHOFF
#define MAT_NEOHOOKEAN
#define MAT_HOOKE
static MoFEMErrorCode setBlocks(MoFEM::Interface &m_field, boost::shared_ptr< map< int, BlockData > > &block_sets_ptr)
Implementation of linear elastic material.
#define MOFEM_LOG_C(channel, severity, format,...)
Implementation of Neo-Hookean elastic material.
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
#define CHKERRG(n)
Check error code of MoFEM/MOAB/PETSc function.
@ BODYFORCESSET
block name is "BODY_FORCES"
@ MAT_ELASTICSET
block name is "MAT_ELASTIC"
@ BLOCKSET
@ MOFEM_NOT_FOUND
Definition definitions.h:33
@ MOFEM_NOT_IMPLEMENTED
Definition definitions.h:32
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define CHKERR
Inline error check.
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
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 bool check_field(const std::string &name) const =0
check if field is in database
#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.
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Manage setting parameters and constitutive equations for nonlinear/linear elastic materials.
std::map< std::string, boost::shared_ptr< NonlinearElasticElement::FunctionsToCalculatePiolaKirchhoffI< adouble > > > aDoubleMaterialModel
virtual MoFEMErrorCode readConfigFile()
read Elastic materials declaration for blocks and meshsets
std::map< std::string, boost::shared_ptr< NonlinearElasticElement::FunctionsToCalculatePiolaKirchhoffI< double > > > doubleMaterialModel
Hash map of materials for evaluation with double.
virtual MoFEMErrorCode iNit()
MoFEMErrorCode setBlocksOrder()
MoFEM::Interface & mField
ElasticMaterials(MoFEM::Interface &m_field, std::string def_material="KIRCHHOFF")
std::map< int, BlockOptionData > blockData
Material parameters on blocks.
po::variables_map vM
PetscBool isConfigFileSet
True if config file is set from command line.
bool iNitialized
true if class is initialized
virtual moab::Interface & get_moab()=0
virtual MPI_Comm & get_comm() const =0
Deprecated interface functions.
Implementation of elastic (non-linear) St. Kirchhoff equation.