v0.15.0
Loading...
Searching...
No Matches
ElasticityMixedFormulation.hpp
Go to the documentation of this file.
1/**
2 * \file ElasticityMixedFormulation.hpp
3 * \example ElasticityMixedFormulation.hpp
4 *
5 * \brief Operator implementation of U-P (mixed) finite element.
6 *
7 */
8
9#ifndef __ELASTICITYMIXEDFORMULATION_HPP__
10#define __ELASTICITYMIXEDFORMULATION_HPP__
11
12struct BlockData {
13 int iD;
14 int oRder;
15 double yOung;
16 double pOisson;
18 BlockData() : oRder(-1), yOung(-1), pOisson(-2) {}
19};
21
22 boost::shared_ptr<MatrixDouble> gradDispPtr;
23 boost::shared_ptr<VectorDouble> pPtr;
25
26 double pOisson;
27 double yOung;
28 double lAmbda;
29 double mU;
30
31 std::map<int, BlockData> setOfBlocksData;
32
34
35 // Setting default values for coeffcients
36 gradDispPtr = boost::shared_ptr<MatrixDouble>(new MatrixDouble());
37 pPtr = boost::shared_ptr<VectorDouble>(new VectorDouble());
38
39 ierr = setBlocks();
40 CHKERRABORT(PETSC_COMM_WORLD, ierr);
41 }
42
44 MoFEMFunctionBegin; // They will be overwriten by BlockData
45 PetscOptionsBegin(PETSC_COMM_WORLD, "", "Problem", "none");
46
47 PetscOptionsEnd();
49 }
50
53
54 yOung = data.yOung;
55 pOisson = data.pOisson;
56 lAmbda = (yOung * pOisson) / ((1. + pOisson) * (1. - 2. * pOisson));
57 mU = yOung / (2. * (1. + pOisson));
58
59 FTensor::Index<'i', 3> i;
60 FTensor::Index<'j', 3> j;
61 FTensor::Index<'k', 3> k;
62 FTensor::Index<'l', 3> l;
63
64 tD(i, j, k, l) = 0.;
65
66 tD(0, 0, 0, 0) = 2 * mU;
67 tD(0, 1, 0, 1) = mU;
68 tD(0, 1, 1, 0) = mU;
69 tD(0, 2, 0, 2) = mU;
70 tD(0, 2, 2, 0) = mU;
71 tD(1, 0, 0, 1) = mU;
72 tD(1, 0, 1, 0) = mU;
73 tD(1, 1, 1, 1) = 2 * mU;
74 tD(1, 2, 1, 2) = mU;
75 tD(1, 2, 2, 1) = mU;
76 tD(2, 0, 0, 2) = mU;
77 tD(2, 0, 2, 0) = mU;
78 tD(2, 1, 1, 2) = mU;
79 tD(2, 1, 2, 1) = mU;
80 tD(2, 2, 2, 2) = 2 * mU;
81
83 }
84
89 Mat_Elastic mydata;
90 CHKERR it->getAttributeDataStructure(mydata);
91 int id = it->getMeshsetId();
92 EntityHandle meshset = it->getMeshset();
93 CHKERR mField.get_moab().get_entities_by_type(
94 meshset, MBTET, setOfBlocksData[id].tEts, true);
95 setOfBlocksData[id].iD = id;
96 setOfBlocksData[id].yOung = mydata.data.Young;
97 setOfBlocksData[id].pOisson = mydata.data.Poisson;
98 }
100 }
101
102private:
104};
105
106/** * @brief Assemble P *
107 * \f[
108 * {\bf{P}} = - \int\limits_\Omega {{\bf{N}}_p^T\frac{1}{\lambda
109 }{{\bf{N}}_p}d\Omega }
110 * \f]
111 *
112 */
114 : public VolumeElementForcesAndSourcesCore::UserDataOperator {
115
117 MatrixDouble locP;
118 MatrixDouble translocP;
120
122 : VolumeElementForcesAndSourcesCore::UserDataOperator(
123 "P", "P", UserDataOperator::OPROWCOL),
124 commonData(common_data), dAta(data) {
125 sYmm = true;
126 }
127
128 PetscErrorCode doWork(int row_side, int col_side, EntityType row_type,
129 EntityType col_type,
130 EntitiesFieldData::EntData &row_data,
131 EntitiesFieldData::EntData &col_data) {
132
134 const int row_nb_dofs = row_data.getIndices().size();
135 if (!row_nb_dofs)
137 const int col_nb_dofs = col_data.getIndices().size();
138 if (!col_nb_dofs)
140
141 if (dAta.tEts.find(getNumeredEntFiniteElementPtr()->getEnt()) ==
142 dAta.tEts.end()) {
144 }
146 // Set size can clear local tangent matrix
147 locP.resize(row_nb_dofs, col_nb_dofs, false);
148 locP.clear();
149
150 const int row_nb_gauss_pts = row_data.getN().size1();
151 if (!row_nb_gauss_pts)
153 const int row_nb_base_functions = row_data.getN().size2();
154 auto row_base_functions = row_data.getFTensor0N();
155
156 // get data
157 FTensor::Index<'i', 3> i;
158 FTensor::Index<'j', 3> j;
159 const double lambda = commonData.lAmbda;
160
161 double coefficient = commonData.pOisson == 0.5 ? 0. : 1 / lambda;
162
163 // integration
164 if (coefficient != 0.) {
165 for (int gg = 0; gg != row_nb_gauss_pts; gg++) {
166
167 // Get volume and integration weight
168 double w = getVolume() * getGaussPts()(3, gg);
169
170 // INTEGRATION
171 int row_bb = 0;
172 for (; row_bb != row_nb_dofs; row_bb++) {
173 auto col_base_functions = col_data.getFTensor0N(gg, 0);
174 for (int col_bb = 0; col_bb != col_nb_dofs; col_bb++) {
175
176 locP(row_bb, col_bb) -=
177 w * row_base_functions * col_base_functions * coefficient;
178
179 ++col_base_functions;
180 }
181 ++row_base_functions;
182 }
183 for (; row_bb != row_nb_base_functions; row_bb++) {
184 ++row_base_functions;
185 }
186 }
187 }
188
189
190 CHKERR MatSetValues(
191 getFEMethod()->ksp_B, row_nb_dofs, &*row_data.getIndices().begin(),
192 col_nb_dofs, &*col_data.getIndices().begin(), &locP(0, 0), ADD_VALUES);
193
194 // is symmetric
195 if (row_side != col_side || row_type != col_type) {
196 translocP.resize(col_nb_dofs, row_nb_dofs, false);
197 noalias(translocP) = trans(locP);
198 CHKERR MatSetValues(getFEMethod()->ksp_B, col_nb_dofs,
199 &*col_data.getIndices().begin(), row_nb_dofs,
200 &*row_data.getIndices().begin(), &translocP(0, 0),
201 ADD_VALUES);
202 }
203
205 }
206};
207
208/** * @brief Assemble G *
209 * \f[
210 * {\bf{G}} = - \int\limits_\Omega {{{\bf{B}}^T}{\bf m}{{\bf{N}}_p}d\Omega }
211 * \f]
212 *
213 */
215 : public VolumeElementForcesAndSourcesCore::UserDataOperator {
216
218 MatrixDouble locG;
220
222 : VolumeElementForcesAndSourcesCore::UserDataOperator(
223 "U", "P", UserDataOperator::OPROWCOL),
224 commonData(common_data), dAta(data) {
225 sYmm = false;
226 }
227
228 PetscErrorCode doWork(int row_side, int col_side, EntityType row_type,
229 EntityType col_type,
230 EntitiesFieldData::EntData &row_data,
231 EntitiesFieldData::EntData &col_data) {
232
234
235 const int row_nb_dofs = row_data.getIndices().size();
236 if (!row_nb_dofs)
238 const int col_nb_dofs = col_data.getIndices().size();
239 if (!col_nb_dofs)
241
242 if (dAta.tEts.find(getNumeredEntFiniteElementPtr()->getEnt()) ==
243 dAta.tEts.end()) {
245 }
247
248 // Set size can clear local tangent matrix
249 locG.resize(row_nb_dofs, col_nb_dofs, false);
250 locG.clear();
251 const int row_nb_gauss_pts = row_data.getN().size1();
252 if (!row_nb_gauss_pts)
254 const int row_nb_base_functions = row_data.getN().size2();
255 auto row_diff_base_functions = row_data.getFTensor1DiffN<3>();
256
258 FTensor::Index<'i', 3> i;
259 FTensor::Index<'j', 3> j;
260
261 // INTEGRATION
262 for (int gg = 0; gg != row_nb_gauss_pts; gg++) {
263
264 // Get volume and integration weight
265 double w = getVolume() * getGaussPts()(3, gg);
266
267 int row_bb = 0;
268 for (; row_bb != row_nb_dofs / 3; row_bb++) {
269
270 t1(i) = w * row_diff_base_functions(i);
271
272 auto base_functions = col_data.getFTensor0N(gg, 0);
273 for (int col_bb = 0; col_bb != col_nb_dofs; col_bb++) {
274
275 FTensor::Tensor1<double *, 3> k(&locG(3 * row_bb + 0, col_bb),
276 &locG(3 * row_bb + 1, col_bb),
277 &locG(3 * row_bb + 2, col_bb));
278
279 k(i) += t1(i) * base_functions;
280 ++base_functions;
281 }
282 ++row_diff_base_functions;
283 }
284 for (; row_bb != row_nb_base_functions; row_bb++) {
285 ++row_diff_base_functions;
286 }
287 }
288
289 CHKERR MatSetValues(getFEMethod()->ksp_B, row_nb_dofs,
290 &*row_data.getIndices().begin(), col_nb_dofs,
291 &*col_data.getIndices().begin(), &*locG.data().begin(),
292 ADD_VALUES);
293
294 // ASSEMBLE THE TRANSPOSE
295 locG = trans(locG);
296 CHKERR MatSetValues(getFEMethod()->ksp_B, col_nb_dofs,
297 &*col_data.getIndices().begin(), row_nb_dofs,
298 &*row_data.getIndices().begin(), &*locG.data().begin(),
299 ADD_VALUES);
301 }
302};
303
304/** * @brief Assemble K *
305 * \f[
306 * {\bf{K}} = \int\limits_\Omega {{{\bf{B}}^T}{{\bf{D}}_d}{\bf{B}}d\Omega }
307 * \f]
308 *
309 */
311 : public VolumeElementForcesAndSourcesCore::UserDataOperator {
312
313 MatrixDouble locK;
314 MatrixDouble translocK;
316
319
321 bool symm = true)
322 : VolumeElementForcesAndSourcesCore::UserDataOperator("U", "U", OPROWCOL,
323 symm),
324 commonData(common_data), dAta(data) {}
325
326 PetscErrorCode doWork(int row_side, int col_side, EntityType row_type,
327 EntityType col_type,
328 EntitiesFieldData::EntData &row_data,
329 EntitiesFieldData::EntData &col_data) {
331
332 auto get_tensor2 = [](MatrixDouble &m, const int r, const int c) {
334 &m(r + 0, c + 0), &m(r + 0, c + 1), &m(r + 0, c + 2),
335 &m(r + 1, c + 0), &m(r + 1, c + 1), &m(r + 1, c + 2),
336 &m(r + 2, c + 0), &m(r + 2, c + 1), &m(r + 2, c + 2));
337 };
338
339 const int row_nb_dofs = row_data.getIndices().size();
340 if (!row_nb_dofs)
342 const int col_nb_dofs = col_data.getIndices().size();
343 if (!col_nb_dofs)
345 if (dAta.tEts.find(getNumeredEntFiniteElementPtr()->getEnt()) ==
346 dAta.tEts.end()) {
348 }
350
351 const bool diagonal_block =
352 (row_type == col_type) && (row_side == col_side);
353 // get number of integration points
354 // Set size can clear local tangent matrix
355 locK.resize(row_nb_dofs, col_nb_dofs, false);
356 locK.clear();
357
358 const int row_nb_gauss_pts = row_data.getN().size1();
359 const int row_nb_base_functions = row_data.getN().size2();
360
361 auto row_diff_base_functions = row_data.getFTensor1DiffN<3>();
362
363 FTensor::Index<'i', 3> i;
364 FTensor::Index<'j', 3> j;
365 FTensor::Index<'k', 3> k;
366 FTensor::Index<'l', 3> l;
367
368 // integrate local matrix for entity block
369 for (int gg = 0; gg != row_nb_gauss_pts; gg++) {
370
371 // Get volume and integration weight
372 double w = getVolume() * getGaussPts()(3, gg);
373
374 int row_bb = 0;
375 for (; row_bb != row_nb_dofs / 3; row_bb++) {
376
377 auto col_diff_base_functions = col_data.getFTensor1DiffN<3>(gg, 0);
378 const int final_bb = diagonal_block ? row_bb + 1 : col_nb_dofs / 3;
379 int col_bb = 0;
380 for (; col_bb != final_bb; col_bb++) {
381
382 auto t_assemble = get_tensor2(locK, 3 * row_bb, 3 * col_bb);
383
384 diffDiff(j, l) =
385 w * row_diff_base_functions(j) * col_diff_base_functions(l);
386
387 t_assemble(i, k) += diffDiff(j, l) * commonData.tD(i, j, k, l);
388 // Next base function for column
389 ++col_diff_base_functions;
390 }
391
392 ++row_diff_base_functions;
393 }
394 for (; row_bb != row_nb_base_functions; row_bb++) {
395 ++row_diff_base_functions;
396 }
397 }
398
399 if (diagonal_block) {
400 for (int row_bb = 0; row_bb != row_nb_dofs / 3; row_bb++) {
401 int col_bb = 0;
402 for (; col_bb != row_bb + 1; col_bb++) {
403 auto t_assemble = get_tensor2(locK, 3 * row_bb, 3 * col_bb);
404 auto t_off_side = get_tensor2(locK, 3 * col_bb, 3 * row_bb);
405 t_off_side(i, k) = t_assemble(k, i);
406 }
407 }
408 }
409
410 const int *row_ind = &*row_data.getIndices().begin();
411 const int *col_ind = &*col_data.getIndices().begin();
412 Mat B = getFEMethod()->ksp_B != PETSC_NULLPTR ? getFEMethod()->ksp_B
413 : getFEMethod()->ksp_B;
414 CHKERR MatSetValues(B, row_nb_dofs, row_ind, col_nb_dofs, col_ind,
415 &locK(0, 0), ADD_VALUES);
416
417 if (row_type != col_type || row_side != col_side) {
418 translocK.resize(col_nb_dofs, row_nb_dofs, false);
419 noalias(translocK) = trans(locK);
420 CHKERR MatSetValues(B, col_nb_dofs, col_ind, row_nb_dofs, row_ind,
421 &translocK(0, 0), ADD_VALUES);
422 }
423
425 }
426};
427
431 moab::Interface &postProcMesh;
432 std::vector<EntityHandle> &mapGaussPts;
434
435 OpPostProcStress(moab::Interface &post_proc_mesh,
436 std::vector<EntityHandle> &map_gauss_pts,
437 DataAtIntegrationPts &common_data, BlockData &data)
438 : VolumeElementForcesAndSourcesCore::UserDataOperator(
439 "U", UserDataOperator::OPROW),
440 commonData(common_data), postProcMesh(post_proc_mesh),
441 mapGaussPts(map_gauss_pts), dAta(data) {
442 doVertices = true;
443 doEdges = false;
444 doQuads = false;
445 doTris = false;
446 doTets = false;
447 doPrisms = false;
448 }
449
450 PetscErrorCode doWork(int side, EntityType type,
451 EntitiesFieldData::EntData &data) {
453 if (type != MBVERTEX)
454 PetscFunctionReturn(9);
455 double def_VAL[9];
456 bzero(def_VAL, 9 * sizeof(double));
457
458 if (dAta.tEts.find(getNumeredEntFiniteElementPtr()->getEnt()) ==
459 dAta.tEts.end()) {
461 }
463
464 Tag th_stress;
465 CHKERR postProcMesh.tag_get_handle("STRESS", 9, MB_TYPE_DOUBLE, th_stress,
466 MB_TAG_CREAT | MB_TAG_SPARSE, def_VAL);
467 Tag th_strain;
468 CHKERR postProcMesh.tag_get_handle("STRAIN", 9, MB_TYPE_DOUBLE, th_strain,
469 MB_TAG_CREAT | MB_TAG_SPARSE, def_VAL);
470 Tag th_psi;
471 CHKERR postProcMesh.tag_get_handle("ENERGY", 1, MB_TYPE_DOUBLE, th_psi,
472 MB_TAG_CREAT | MB_TAG_SPARSE, def_VAL);
473
474 auto grad = getFTensor2FromMat<3, 3>(*commonData.gradDispPtr);
475 auto p = getFTensor0FromVec(*commonData.pPtr);
476
477 const int nb_gauss_pts = commonData.gradDispPtr->size2();
478 const double mu = commonData.mU;
479
480 FTensor::Index<'i', 3> i;
481 FTensor::Index<'j', 3> j;
484
485 for (int gg = 0; gg != nb_gauss_pts; gg++) {
486 strain(i, j) = 0.5 * (grad(i, j) + grad(j, i));
487 double psi = 0.5 * p * p + mu * strain(i, j) * strain(i, j);
488
489 stress(i, j) = 2 * mu * strain(i, j);
490 stress(1, 1) -= p;
491 stress(0, 0) -= p;
492 stress(2, 2) -= p;
493
494 CHKERR postProcMesh.tag_set_data(th_psi, &mapGaussPts[gg], 1, &psi);
495 CHKERR postProcMesh.tag_set_data(th_strain, &mapGaussPts[gg], 1,
496 &strain(0, 0));
497 CHKERR postProcMesh.tag_set_data(th_stress, &mapGaussPts[gg], 1,
498 &stress(0, 0));
499 ++p;
500 ++grad;
501 }
502
504 }
505};
506
507#endif //__ELASTICITYMIXEDFORMULATION_HPP__
ForcesAndSourcesCore::UserDataOperator UserDataOperator
#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 ...
@ MAT_ELASTICSET
block name is "MAT_ELASTIC"
@ BLOCKSET
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define CHKERR
Inline error check.
#define _IT_CUBITMESHSETS_BY_BCDATA_TYPE_FOR_LOOP_(MESHSET_MANAGER, CUBITBCTYPE, IT)
Iterator that loops over a specific Cubit MeshSet in a moFEM field.
FTensor::Index< 'i', SPACE_DIM > i
static double lambda
const double c
speed of light (cm/ns)
FTensor::Index< 'l', 3 > l
FTensor::Index< 'j', 3 > j
FTensor::Index< 'k', 3 > k
static MoFEMErrorCodeGeneric< PetscErrorCode > ierr
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
UBlasMatrix< double > MatrixDouble
Definition Types.hpp:77
UBlasVector< double > VectorDouble
Definition Types.hpp:68
FTensor::Index< 'm', 3 > m
MoFEMErrorCode getBlockData(BlockData &data)
boost::shared_ptr< MatrixDouble > gradDispPtr
FTensor::Ddg< double, 3, 3 > tD
std::map< int, BlockData > setOfBlocksData
DataAtIntegrationPts(MoFEM::Interface &m_field)
boost::shared_ptr< VectorDouble > pPtr
virtual moab::Interface & get_moab()=0
bool & doTris
\deprectaed
bool & doPrisms
\deprectaed
bool & doVertices
\deprectaed If false skip vertices
bool & doTets
\deprectaed
bool & doEdges
\deprectaed If false skip edges
bool & doQuads
\deprectaed
Deprecated interface functions.
boost::shared_ptr< const NumeredEntFiniteElement > getNumeredEntFiniteElementPtr() const
Return raw pointer to NumeredEntFiniteElement.
@ OPROW
operator doWork function is executed on FE rows
Elastic material data structure.
OpAssembleG(DataAtIntegrationPts &common_data, BlockData &data)
PetscErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
DataAtIntegrationPts & commonData
OpAssembleK(DataAtIntegrationPts &common_data, BlockData &data, bool symm=true)
DataAtIntegrationPts & commonData
PetscErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
FTensor::Tensor2< double, 3, 3 > diffDiff
DataAtIntegrationPts & commonData
PetscErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
OpAssembleP(DataAtIntegrationPts &common_data, BlockData &data)
std::vector< EntityHandle > & mapGaussPts
OpPostProcStress(moab::Interface &post_proc_mesh, std::vector< EntityHandle > &map_gauss_pts, DataAtIntegrationPts &common_data, BlockData &data)
DataAtIntegrationPts & commonData
PetscErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)