v0.14.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
43 MoFEMErrorCode getParameters() {
44 MoFEMFunctionBegin; // They will be overwriten by BlockData
45 CHKERR PetscOptionsBegin(PETSC_COMM_WORLD, "", "Problem", "none");
46
47 ierr = PetscOptionsEnd();
48 CHKERRQ(ierr);
50 }
51
52 MoFEMErrorCode getBlockData(BlockData &data) {
54
55 yOung = data.yOung;
56 pOisson = data.pOisson;
57 lAmbda = (yOung * pOisson) / ((1. + pOisson) * (1. - 2. * pOisson));
58 mU = yOung / (2. * (1. + pOisson));
59
64
65 tD(i, j, k, l) = 0.;
66
67 tD(0, 0, 0, 0) = 2 * mU;
68 tD(0, 1, 0, 1) = mU;
69 tD(0, 1, 1, 0) = mU;
70 tD(0, 2, 0, 2) = mU;
71 tD(0, 2, 2, 0) = mU;
72 tD(1, 0, 0, 1) = mU;
73 tD(1, 0, 1, 0) = mU;
74 tD(1, 1, 1, 1) = 2 * mU;
75 tD(1, 2, 1, 2) = mU;
76 tD(1, 2, 2, 1) = mU;
77 tD(2, 0, 0, 2) = mU;
78 tD(2, 0, 2, 0) = mU;
79 tD(2, 1, 1, 2) = mU;
80 tD(2, 1, 2, 1) = mU;
81 tD(2, 2, 2, 2) = 2 * mU;
82
84 }
85
86 MoFEMErrorCode setBlocks() {
90 Mat_Elastic mydata;
91 CHKERR it->getAttributeDataStructure(mydata);
92 int id = it->getMeshsetId();
93 EntityHandle meshset = it->getMeshset();
94 CHKERR mField.get_moab().get_entities_by_type(
95 meshset, MBTET, setOfBlocksData[id].tEts, true);
96 setOfBlocksData[id].iD = id;
97 setOfBlocksData[id].yOung = mydata.data.Young;
98 setOfBlocksData[id].pOisson = mydata.data.Poisson;
99 }
101 }
102
103private:
105};
106
107/** * @brief Assemble P *
108 * \f[
109 * {\bf{P}} = - \int\limits_\Omega {{\bf{N}}_p^T\frac{1}{\lambda
110 }{{\bf{N}}_p}d\Omega }
111 * \f]
112 *
113 */
115 : public VolumeElementForcesAndSourcesCore::UserDataOperator {
116
118 MatrixDouble locP;
119 MatrixDouble translocP;
121
123 : VolumeElementForcesAndSourcesCore::UserDataOperator(
124 "P", "P", UserDataOperator::OPROWCOL),
125 commonData(common_data), dAta(data) {
126 sYmm = true;
127 }
128
129 PetscErrorCode doWork(int row_side, int col_side, EntityType row_type,
130 EntityType col_type,
131 EntitiesFieldData::EntData &row_data,
132 EntitiesFieldData::EntData &col_data) {
133
135 const int row_nb_dofs = row_data.getIndices().size();
136 if (!row_nb_dofs)
138 const int col_nb_dofs = col_data.getIndices().size();
139 if (!col_nb_dofs)
141
142 if (dAta.tEts.find(getNumeredEntFiniteElementPtr()->getEnt()) ==
143 dAta.tEts.end()) {
145 }
147 // Set size can clear local tangent matrix
148 locP.resize(row_nb_dofs, col_nb_dofs, false);
149 locP.clear();
150
151 const int row_nb_gauss_pts = row_data.getN().size1();
152 if (!row_nb_gauss_pts)
154 const int row_nb_base_functions = row_data.getN().size2();
155 auto row_base_functions = row_data.getFTensor0N();
156
157 // get data
160 const double lambda = commonData.lAmbda;
161
162 double coefficient = commonData.pOisson == 0.5 ? 0. : 1 / lambda;
163
164 // integration
165 if (coefficient != 0.) {
166 for (int gg = 0; gg != row_nb_gauss_pts; gg++) {
167
168 // Get volume and integration weight
169 double w = getVolume() * getGaussPts()(3, gg);
170
171 // INTEGRATION
172 int row_bb = 0;
173 for (; row_bb != row_nb_dofs; row_bb++) {
174 auto col_base_functions = col_data.getFTensor0N(gg, 0);
175 for (int col_bb = 0; col_bb != col_nb_dofs; col_bb++) {
176
177 locP(row_bb, col_bb) -=
178 w * row_base_functions * col_base_functions * coefficient;
179
180 ++col_base_functions;
181 }
182 ++row_base_functions;
183 }
184 for (; row_bb != row_nb_base_functions; row_bb++) {
185 ++row_base_functions;
186 }
187 }
188 }
189
190
191 CHKERR MatSetValues(
192 getFEMethod()->ksp_B, row_nb_dofs, &*row_data.getIndices().begin(),
193 col_nb_dofs, &*col_data.getIndices().begin(), &locP(0, 0), ADD_VALUES);
194
195 // is symmetric
196 if (row_side != col_side || row_type != col_type) {
197 translocP.resize(col_nb_dofs, row_nb_dofs, false);
198 noalias(translocP) = trans(locP);
199 CHKERR MatSetValues(getFEMethod()->ksp_B, col_nb_dofs,
200 &*col_data.getIndices().begin(), row_nb_dofs,
201 &*row_data.getIndices().begin(), &translocP(0, 0),
202 ADD_VALUES);
203 }
204
206 }
207};
208
209/** * @brief Assemble G *
210 * \f[
211 * {\bf{G}} = - \int\limits_\Omega {{{\bf{B}}^T}{\bf m}{{\bf{N}}_p}d\Omega }
212 * \f]
213 *
214 */
216 : public VolumeElementForcesAndSourcesCore::UserDataOperator {
217
219 MatrixDouble locG;
221
223 : VolumeElementForcesAndSourcesCore::UserDataOperator(
224 "U", "P", UserDataOperator::OPROWCOL),
225 commonData(common_data), dAta(data) {
226 sYmm = false;
227 }
228
229 PetscErrorCode doWork(int row_side, int col_side, EntityType row_type,
230 EntityType col_type,
231 EntitiesFieldData::EntData &row_data,
232 EntitiesFieldData::EntData &col_data) {
233
235
236 const int row_nb_dofs = row_data.getIndices().size();
237 if (!row_nb_dofs)
239 const int col_nb_dofs = col_data.getIndices().size();
240 if (!col_nb_dofs)
242
243 if (dAta.tEts.find(getNumeredEntFiniteElementPtr()->getEnt()) ==
244 dAta.tEts.end()) {
246 }
248
249 // Set size can clear local tangent matrix
250 locG.resize(row_nb_dofs, col_nb_dofs, false);
251 locG.clear();
252 const int row_nb_gauss_pts = row_data.getN().size1();
253 if (!row_nb_gauss_pts)
255 const int row_nb_base_functions = row_data.getN().size2();
256 auto row_diff_base_functions = row_data.getFTensor1DiffN<3>();
257
261
262 // INTEGRATION
263 for (int gg = 0; gg != row_nb_gauss_pts; gg++) {
264
265 // Get volume and integration weight
266 double w = getVolume() * getGaussPts()(3, gg);
267
268 int row_bb = 0;
269 for (; row_bb != row_nb_dofs / 3; row_bb++) {
270
271 t1(i) = w * row_diff_base_functions(i);
272
273 auto base_functions = col_data.getFTensor0N(gg, 0);
274 for (int col_bb = 0; col_bb != col_nb_dofs; col_bb++) {
275
276 FTensor::Tensor1<double *, 3> k(&locG(3 * row_bb + 0, col_bb),
277 &locG(3 * row_bb + 1, col_bb),
278 &locG(3 * row_bb + 2, col_bb));
279
280 k(i) += t1(i) * base_functions;
281 ++base_functions;
282 }
283 ++row_diff_base_functions;
284 }
285 for (; row_bb != row_nb_base_functions; row_bb++) {
286 ++row_diff_base_functions;
287 }
288 }
289
290 CHKERR MatSetValues(getFEMethod()->ksp_B, row_nb_dofs,
291 &*row_data.getIndices().begin(), col_nb_dofs,
292 &*col_data.getIndices().begin(), &*locG.data().begin(),
293 ADD_VALUES);
294
295 // ASSEMBLE THE TRANSPOSE
296 locG = trans(locG);
297 CHKERR MatSetValues(getFEMethod()->ksp_B, col_nb_dofs,
298 &*col_data.getIndices().begin(), row_nb_dofs,
299 &*row_data.getIndices().begin(), &*locG.data().begin(),
300 ADD_VALUES);
302 }
303};
304
305/** * @brief Assemble K *
306 * \f[
307 * {\bf{K}} = \int\limits_\Omega {{{\bf{B}}^T}{{\bf{D}}_d}{\bf{B}}d\Omega }
308 * \f]
309 *
310 */
311struct OpAssembleK
312 : public VolumeElementForcesAndSourcesCore::UserDataOperator {
313
314 MatrixDouble locK;
315 MatrixDouble translocK;
317
320
322 bool symm = true)
323 : VolumeElementForcesAndSourcesCore::UserDataOperator("U", "U", OPROWCOL,
324 symm),
325 commonData(common_data), dAta(data) {}
326
327 PetscErrorCode doWork(int row_side, int col_side, EntityType row_type,
328 EntityType col_type,
329 EntitiesFieldData::EntData &row_data,
330 EntitiesFieldData::EntData &col_data) {
332
333 auto get_tensor2 = [](MatrixDouble &m, const int r, const int c) {
335 &m(r + 0, c + 0), &m(r + 0, c + 1), &m(r + 0, c + 2),
336 &m(r + 1, c + 0), &m(r + 1, c + 1), &m(r + 1, c + 2),
337 &m(r + 2, c + 0), &m(r + 2, c + 1), &m(r + 2, c + 2));
338 };
339
340 const int row_nb_dofs = row_data.getIndices().size();
341 if (!row_nb_dofs)
343 const int col_nb_dofs = col_data.getIndices().size();
344 if (!col_nb_dofs)
346 if (dAta.tEts.find(getNumeredEntFiniteElementPtr()->getEnt()) ==
347 dAta.tEts.end()) {
349 }
350 commonData.getBlockData(dAta);
351
352 const bool diagonal_block =
353 (row_type == col_type) && (row_side == col_side);
354 // get number of integration points
355 // Set size can clear local tangent matrix
356 locK.resize(row_nb_dofs, col_nb_dofs, false);
357 locK.clear();
358
359 const int row_nb_gauss_pts = row_data.getN().size1();
360 const int row_nb_base_functions = row_data.getN().size2();
361
362 auto row_diff_base_functions = row_data.getFTensor1DiffN<3>();
363
368
369 // integrate local matrix for entity block
370 for (int gg = 0; gg != row_nb_gauss_pts; gg++) {
371
372 // Get volume and integration weight
373 double w = getVolume() * getGaussPts()(3, gg);
374
375 int row_bb = 0;
376 for (; row_bb != row_nb_dofs / 3; row_bb++) {
377
378 auto col_diff_base_functions = col_data.getFTensor1DiffN<3>(gg, 0);
379 const int final_bb = diagonal_block ? row_bb + 1 : col_nb_dofs / 3;
380 int col_bb = 0;
381 for (; col_bb != final_bb; col_bb++) {
382
383 auto t_assemble = get_tensor2(locK, 3 * row_bb, 3 * col_bb);
384
385 diffDiff(j, l) =
386 w * row_diff_base_functions(j) * col_diff_base_functions(l);
387
388 t_assemble(i, k) += diffDiff(j, l) * commonData.tD(i, j, k, l);
389 // Next base function for column
390 ++col_diff_base_functions;
391 }
392
393 ++row_diff_base_functions;
394 }
395 for (; row_bb != row_nb_base_functions; row_bb++) {
396 ++row_diff_base_functions;
397 }
398 }
399
400 if (diagonal_block) {
401 for (int row_bb = 0; row_bb != row_nb_dofs / 3; row_bb++) {
402 int col_bb = 0;
403 for (; col_bb != row_bb + 1; col_bb++) {
404 auto t_assemble = get_tensor2(locK, 3 * row_bb, 3 * col_bb);
405 auto t_off_side = get_tensor2(locK, 3 * col_bb, 3 * row_bb);
406 t_off_side(i, k) = t_assemble(k, i);
407 }
408 }
409 }
410
411 const int *row_ind = &*row_data.getIndices().begin();
412 const int *col_ind = &*col_data.getIndices().begin();
413 Mat B = getFEMethod()->ksp_B != PETSC_NULL ? getFEMethod()->ksp_B
414 : getFEMethod()->ksp_B;
415 CHKERR MatSetValues(B, row_nb_dofs, row_ind, col_nb_dofs, col_ind,
416 &locK(0, 0), ADD_VALUES);
417
418 if (row_type != col_type || row_side != col_side) {
419 translocK.resize(col_nb_dofs, row_nb_dofs, false);
420 noalias(translocK) = trans(locK);
421 CHKERR MatSetValues(B, col_nb_dofs, col_ind, row_nb_dofs, row_ind,
422 &translocK(0, 0), ADD_VALUES);
423 }
424
426 }
427};
428
432 moab::Interface &postProcMesh;
433 std::vector<EntityHandle> &mapGaussPts;
435
436 OpPostProcStress(moab::Interface &post_proc_mesh,
437 std::vector<EntityHandle> &map_gauss_pts,
438 DataAtIntegrationPts &common_data, BlockData &data)
439 : VolumeElementForcesAndSourcesCore::UserDataOperator(
440 "U", UserDataOperator::OPROW),
441 commonData(common_data), postProcMesh(post_proc_mesh),
442 mapGaussPts(map_gauss_pts), dAta(data) {
443 doVertices = true;
444 doEdges = false;
445 doQuads = false;
446 doTris = false;
447 doTets = false;
448 doPrisms = false;
449 }
450
451 PetscErrorCode doWork(int side, EntityType type,
452 EntitiesFieldData::EntData &data) {
454 if (type != MBVERTEX)
455 PetscFunctionReturn(9);
456 double def_VAL[9];
457 bzero(def_VAL, 9 * sizeof(double));
458
459 if (dAta.tEts.find(getNumeredEntFiniteElementPtr()->getEnt()) ==
460 dAta.tEts.end()) {
462 }
464
465 Tag th_stress;
466 CHKERR postProcMesh.tag_get_handle("STRESS", 9, MB_TYPE_DOUBLE, th_stress,
467 MB_TAG_CREAT | MB_TAG_SPARSE, def_VAL);
468 Tag th_strain;
469 CHKERR postProcMesh.tag_get_handle("STRAIN", 9, MB_TYPE_DOUBLE, th_strain,
470 MB_TAG_CREAT | MB_TAG_SPARSE, def_VAL);
471 Tag th_psi;
472 CHKERR postProcMesh.tag_get_handle("ENERGY", 1, MB_TYPE_DOUBLE, th_psi,
473 MB_TAG_CREAT | MB_TAG_SPARSE, def_VAL);
474
475 auto grad = getFTensor2FromMat<3, 3>(*commonData.gradDispPtr);
476 auto p = getFTensor0FromVec(*commonData.pPtr);
477
478 const int nb_gauss_pts = commonData.gradDispPtr->size2();
479 const double mu = commonData.mU;
480
485
486 for (int gg = 0; gg != nb_gauss_pts; gg++) {
487 strain(i, j) = 0.5 * (grad(i, j) + grad(j, i));
488 double psi = 0.5 * p * p + mu * strain(i, j) * strain(i, j);
489
490 stress(i, j) = 2 * mu * strain(i, j);
491 stress(1, 1) -= p;
492 stress(0, 0) -= p;
493 stress(2, 2) -= p;
494
495 CHKERR postProcMesh.tag_set_data(th_psi, &mapGaussPts[gg], 1, &psi);
496 CHKERR postProcMesh.tag_set_data(th_strain, &mapGaussPts[gg], 1,
497 &strain(0, 0));
498 CHKERR postProcMesh.tag_set_data(th_stress, &mapGaussPts[gg], 1,
499 &stress(0, 0));
500 ++p;
501 ++grad;
502 }
503
505 }
506};
507
508#endif //__ELASTICITYMIXEDFORMULATION_HPP__
static Index< 'p', 3 > p
ForcesAndSourcesCore::UserDataOperator UserDataOperator
static PetscErrorCode ierr
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:447
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:346
@ MAT_ELASTICSET
block name is "MAT_ELASTIC"
Definition: definitions.h:159
@ BLOCKSET
Definition: definitions.h:148
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:416
#define CHKERR
Inline error check.
Definition: definitions.h:535
FTensor::Index< 'm', SPACE_DIM > m
#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
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
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
Assemble K *.
Definition: Header.hpp:339
CommonData & commonData
Definition: Header.hpp:346
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
BlockData & dAta
moab::Interface & postProcMesh
PetscErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)