9#ifndef __ELASTICITYMIXEDFORMULATION_HPP__
10#define __ELASTICITYMIXEDFORMULATION_HPP__
23 boost::shared_ptr<VectorDouble>
pPtr;
36 gradDispPtr = boost::shared_ptr<MatrixDouble>(
new MatrixDouble());
37 pPtr = boost::shared_ptr<VectorDouble>(
new VectorDouble());
40 CHKERRABORT(PETSC_COMM_WORLD,
ierr);
45 CHKERR PetscOptionsBegin(PETSC_COMM_WORLD,
"",
"Problem",
"none");
47 ierr = PetscOptionsEnd();
67 tD(0, 0, 0, 0) = 2 *
mU;
74 tD(1, 1, 1, 1) = 2 *
mU;
81 tD(2, 2, 2, 2) = 2 *
mU;
91 CHKERR it->getAttributeDataStructure(mydata);
92 int id = it->getMeshsetId();
115 :
public VolumeElementForcesAndSourcesCore::UserDataOperator {
131 EntitiesFieldData::EntData &row_data,
132 EntitiesFieldData::EntData &col_data) {
135 const int row_nb_dofs = row_data.getIndices().size();
138 const int col_nb_dofs = col_data.getIndices().size();
142 if (
dAta.
tEts.find(getNumeredEntFiniteElementPtr()->getEnt()) ==
148 locP.resize(row_nb_dofs, col_nb_dofs,
false);
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();
165 if (coefficient != 0.) {
166 for (
int gg = 0; gg != row_nb_gauss_pts; gg++) {
169 double w = getVolume() * getGaussPts()(3, gg);
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++) {
177 locP(row_bb, col_bb) -=
178 w * row_base_functions * col_base_functions * coefficient;
180 ++col_base_functions;
182 ++row_base_functions;
184 for (; row_bb != row_nb_base_functions; row_bb++) {
185 ++row_base_functions;
192 getFEMethod()->ksp_B, row_nb_dofs, &*row_data.getIndices().begin(),
193 col_nb_dofs, &*col_data.getIndices().begin(), &
locP(0, 0), ADD_VALUES);
196 if (row_side != col_side || row_type != col_type) {
197 translocP.resize(col_nb_dofs, row_nb_dofs,
false);
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),
216 :
public VolumeElementForcesAndSourcesCore::UserDataOperator {
231 EntitiesFieldData::EntData &row_data,
232 EntitiesFieldData::EntData &col_data) {
236 const int row_nb_dofs = row_data.getIndices().size();
239 const int col_nb_dofs = col_data.getIndices().size();
243 if (
dAta.
tEts.find(getNumeredEntFiniteElementPtr()->getEnt()) ==
250 locG.resize(row_nb_dofs, col_nb_dofs,
false);
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>();
263 for (
int gg = 0; gg != row_nb_gauss_pts; gg++) {
266 double w = getVolume() * getGaussPts()(3, gg);
269 for (; row_bb != row_nb_dofs / 3; row_bb++) {
271 t1(
i) = w * row_diff_base_functions(
i);
273 auto base_functions = col_data.getFTensor0N(gg, 0);
274 for (
int col_bb = 0; col_bb != col_nb_dofs; col_bb++) {
277 &
locG(3 * row_bb + 1, col_bb),
278 &
locG(3 * row_bb + 2, col_bb));
280 k(
i) += t1(
i) * base_functions;
283 ++row_diff_base_functions;
285 for (; row_bb != row_nb_base_functions; row_bb++) {
286 ++row_diff_base_functions;
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(),
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(),
312 :
public VolumeElementForcesAndSourcesCore::UserDataOperator {
329 EntitiesFieldData::EntData &row_data,
330 EntitiesFieldData::EntData &col_data) {
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));
340 const int row_nb_dofs = row_data.getIndices().size();
343 const int col_nb_dofs = col_data.getIndices().size();
346 if (
dAta.
tEts.find(getNumeredEntFiniteElementPtr()->getEnt()) ==
352 const bool diagonal_block =
353 (row_type == col_type) && (row_side == col_side);
356 locK.resize(row_nb_dofs, col_nb_dofs,
false);
359 const int row_nb_gauss_pts = row_data.getN().size1();
360 const int row_nb_base_functions = row_data.getN().size2();
362 auto row_diff_base_functions = row_data.getFTensor1DiffN<3>();
370 for (
int gg = 0; gg != row_nb_gauss_pts; gg++) {
373 double w = getVolume() * getGaussPts()(3, gg);
376 for (; row_bb != row_nb_dofs / 3; row_bb++) {
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;
381 for (; col_bb != final_bb; col_bb++) {
383 auto t_assemble = get_tensor2(
locK, 3 * row_bb, 3 * col_bb);
386 w * row_diff_base_functions(
j) * col_diff_base_functions(
l);
390 ++col_diff_base_functions;
393 ++row_diff_base_functions;
395 for (; row_bb != row_nb_base_functions; row_bb++) {
396 ++row_diff_base_functions;
400 if (diagonal_block) {
401 for (
int row_bb = 0; row_bb != row_nb_dofs / 3; row_bb++) {
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);
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);
418 if (row_type != col_type || row_side != col_side) {
419 translocK.resize(col_nb_dofs, row_nb_dofs,
false);
421 CHKERR MatSetValues(
B, col_nb_dofs, col_ind, row_nb_dofs, row_ind,
437 std::vector<EntityHandle> &map_gauss_pts,
452 EntitiesFieldData::EntData &data) {
454 if (type != MBVERTEX)
455 PetscFunctionReturn(9);
457 bzero(def_VAL, 9 *
sizeof(
double));
467 MB_TAG_CREAT | MB_TAG_SPARSE, def_VAL);
470 MB_TAG_CREAT | MB_TAG_SPARSE, def_VAL);
473 MB_TAG_CREAT | MB_TAG_SPARSE, def_VAL);
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);
490 stress(
i,
j) = 2 *
mu * strain(
i,
j);
ForcesAndSourcesCore::UserDataOperator UserDataOperator
static PetscErrorCode ierr
#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"
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define CHKERR
Inline error check.
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
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
MoFEM::Interface & mField
std::map< int, BlockData > setOfBlocksData
DataAtIntegrationPts(MoFEM::Interface &m_field)
MoFEMErrorCode setBlocks()
boost::shared_ptr< VectorDouble > pPtr
MoFEMErrorCode getParameters()
virtual moab::Interface & get_moab()=0
bool & doPrisms
\deprectaed
bool & doVertices
\deprectaed If false skip vertices
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
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
moab::Interface & postProcMesh
PetscErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)