v0.14.0
ElasticityMixedFormulation.hpp

Operator implementation of U-P (mixed) finite element.

/**
* \file ElasticityMixedFormulation.hpp
* \example ElasticityMixedFormulation.hpp
*
* \brief Operator implementation of U-P (mixed) finite element.
*
*/
#ifndef __ELASTICITYMIXEDFORMULATION_HPP__
#define __ELASTICITYMIXEDFORMULATION_HPP__
struct BlockData {
int iD;
int oRder;
double yOung;
double pOisson;
BlockData() : oRder(-1), yOung(-1), pOisson(-2) {}
};
boost::shared_ptr<VectorDouble> pPtr;
double pOisson;
double yOung;
double lAmbda;
double mU;
std::map<int, BlockData> setOfBlocksData;
// Setting default values for coeffcients
pPtr = boost::shared_ptr<VectorDouble>(new VectorDouble());
CHKERRABORT(PETSC_COMM_WORLD, ierr);
}
MoFEMFunctionBegin; // They will be overwriten by BlockData
CHKERR PetscOptionsBegin(PETSC_COMM_WORLD, "", "Problem", "none");
ierr = PetscOptionsEnd();
CHKERRQ(ierr);
}
yOung = data.yOung;
pOisson = data.pOisson;
lAmbda = (yOung * pOisson) / ((1. + pOisson) * (1. - 2. * pOisson));
mU = yOung / (2. * (1. + pOisson));
tD(i, j, k, l) = 0.;
tD(0, 0, 0, 0) = 2 * mU;
tD(0, 1, 0, 1) = mU;
tD(0, 1, 1, 0) = mU;
tD(0, 2, 0, 2) = mU;
tD(0, 2, 2, 0) = mU;
tD(1, 0, 0, 1) = mU;
tD(1, 0, 1, 0) = mU;
tD(1, 1, 1, 1) = 2 * mU;
tD(1, 2, 1, 2) = mU;
tD(1, 2, 2, 1) = mU;
tD(2, 0, 0, 2) = mU;
tD(2, 0, 2, 0) = mU;
tD(2, 1, 1, 2) = mU;
tD(2, 1, 2, 1) = mU;
tD(2, 2, 2, 2) = 2 * mU;
}
Mat_Elastic mydata;
CHKERR it->getAttributeDataStructure(mydata);
int id = it->getMeshsetId();
EntityHandle meshset = it->getMeshset();
CHKERR mField.get_moab().get_entities_by_type(
meshset, MBTET, setOfBlocksData[id].tEts, true);
setOfBlocksData[id].iD = id;
setOfBlocksData[id].yOung = mydata.data.Young;
setOfBlocksData[id].pOisson = mydata.data.Poisson;
}
}
private:
};
/** * @brief Assemble P *
* \f[
* {\bf{P}} = - \int\limits_\Omega {{\bf{N}}_p^T\frac{1}{\lambda
}{{\bf{N}}_p}d\Omega }
* \f]
*
*/
: VolumeElementForcesAndSourcesCore::UserDataOperator(
"P", "P", UserDataOperator::OPROWCOL),
commonData(common_data), dAta(data) {
sYmm = true;
}
PetscErrorCode doWork(int row_side, int col_side, EntityType row_type,
EntityType col_type,
const int row_nb_dofs = row_data.getIndices().size();
if (!row_nb_dofs)
const int col_nb_dofs = col_data.getIndices().size();
if (!col_nb_dofs)
if (dAta.tEts.find(getNumeredEntFiniteElementPtr()->getEnt()) ==
dAta.tEts.end()) {
}
// Set size can clear local tangent matrix
locP.resize(row_nb_dofs, col_nb_dofs, false);
locP.clear();
const int row_nb_gauss_pts = row_data.getN().size1();
if (!row_nb_gauss_pts)
const int row_nb_base_functions = row_data.getN().size2();
auto row_base_functions = row_data.getFTensor0N();
// get data
const double lambda = commonData.lAmbda;
double coefficient = commonData.pOisson == 0.5 ? 0. : 1 / lambda;
// integration
if (coefficient != 0.) {
for (int gg = 0; gg != row_nb_gauss_pts; gg++) {
// Get volume and integration weight
double w = getVolume() * getGaussPts()(3, gg);
// INTEGRATION
int row_bb = 0;
for (; row_bb != row_nb_dofs; row_bb++) {
auto col_base_functions = col_data.getFTensor0N(gg, 0);
for (int col_bb = 0; col_bb != col_nb_dofs; col_bb++) {
locP(row_bb, col_bb) -=
w * row_base_functions * col_base_functions * coefficient;
++col_base_functions;
}
++row_base_functions;
}
for (; row_bb != row_nb_base_functions; row_bb++) {
++row_base_functions;
}
}
}
getFEMethod()->ksp_B, row_nb_dofs, &*row_data.getIndices().begin(),
// is symmetric
if (row_side != col_side || row_type != col_type) {
translocP.resize(col_nb_dofs, row_nb_dofs, false);
noalias(translocP) = trans(locP);
CHKERR MatSetValues(getFEMethod()->ksp_B, col_nb_dofs,
&*col_data.getIndices().begin(), row_nb_dofs,
&*row_data.getIndices().begin(), &translocP(0, 0),
}
}
};
/** * @brief Assemble G *
* \f[
* {\bf{G}} = - \int\limits_\Omega {{{\bf{B}}^T}{\bf m}{{\bf{N}}_p}d\Omega }
* \f]
*
*/
: VolumeElementForcesAndSourcesCore::UserDataOperator(
"U", "P", UserDataOperator::OPROWCOL),
commonData(common_data), dAta(data) {
sYmm = false;
}
PetscErrorCode doWork(int row_side, int col_side, EntityType row_type,
EntityType col_type,
const int row_nb_dofs = row_data.getIndices().size();
if (!row_nb_dofs)
const int col_nb_dofs = col_data.getIndices().size();
if (!col_nb_dofs)
if (dAta.tEts.find(getNumeredEntFiniteElementPtr()->getEnt()) ==
dAta.tEts.end()) {
}
// Set size can clear local tangent matrix
locG.resize(row_nb_dofs, col_nb_dofs, false);
locG.clear();
const int row_nb_gauss_pts = row_data.getN().size1();
if (!row_nb_gauss_pts)
const int row_nb_base_functions = row_data.getN().size2();
auto row_diff_base_functions = row_data.getFTensor1DiffN<3>();
// INTEGRATION
for (int gg = 0; gg != row_nb_gauss_pts; gg++) {
// Get volume and integration weight
double w = getVolume() * getGaussPts()(3, gg);
int row_bb = 0;
for (; row_bb != row_nb_dofs / 3; row_bb++) {
t1(i) = w * row_diff_base_functions(i);
auto base_functions = col_data.getFTensor0N(gg, 0);
for (int col_bb = 0; col_bb != col_nb_dofs; col_bb++) {
FTensor::Tensor1<double *, 3> k(&locG(3 * row_bb + 0, col_bb),
&locG(3 * row_bb + 1, col_bb),
&locG(3 * row_bb + 2, col_bb));
k(i) += t1(i) * base_functions;
++base_functions;
}
++row_diff_base_functions;
}
for (; row_bb != row_nb_base_functions; row_bb++) {
++row_diff_base_functions;
}
}
CHKERR MatSetValues(getFEMethod()->ksp_B, row_nb_dofs,
&*row_data.getIndices().begin(), col_nb_dofs,
&*col_data.getIndices().begin(), &*locG.data().begin(),
// ASSEMBLE THE TRANSPOSE
locG = trans(locG);
CHKERR MatSetValues(getFEMethod()->ksp_B, col_nb_dofs,
&*col_data.getIndices().begin(), row_nb_dofs,
&*row_data.getIndices().begin(), &*locG.data().begin(),
}
};
/** * @brief Assemble K *
* \f[
* {\bf{K}} = \int\limits_\Omega {{{\bf{B}}^T}{{\bf{D}}_d}{\bf{B}}d\Omega }
* \f]
*
*/
bool symm = true)
: VolumeElementForcesAndSourcesCore::UserDataOperator("U", "U", OPROWCOL,
symm),
commonData(common_data), dAta(data) {}
PetscErrorCode doWork(int row_side, int col_side, EntityType row_type,
EntityType col_type,
auto get_tensor2 = [](MatrixDouble &m, const int r, const int c) {
&m(r + 0, c + 0), &m(r + 0, c + 1), &m(r + 0, c + 2),
&m(r + 1, c + 0), &m(r + 1, c + 1), &m(r + 1, c + 2),
&m(r + 2, c + 0), &m(r + 2, c + 1), &m(r + 2, c + 2));
};
const int row_nb_dofs = row_data.getIndices().size();
if (!row_nb_dofs)
const int col_nb_dofs = col_data.getIndices().size();
if (!col_nb_dofs)
if (dAta.tEts.find(getNumeredEntFiniteElementPtr()->getEnt()) ==
dAta.tEts.end()) {
}
const bool diagonal_block =
(row_type == col_type) && (row_side == col_side);
// get number of integration points
// Set size can clear local tangent matrix
locK.resize(row_nb_dofs, col_nb_dofs, false);
locK.clear();
const int row_nb_gauss_pts = row_data.getN().size1();
const int row_nb_base_functions = row_data.getN().size2();
auto row_diff_base_functions = row_data.getFTensor1DiffN<3>();
// integrate local matrix for entity block
for (int gg = 0; gg != row_nb_gauss_pts; gg++) {
// Get volume and integration weight
double w = getVolume() * getGaussPts()(3, gg);
int row_bb = 0;
for (; row_bb != row_nb_dofs / 3; row_bb++) {
auto col_diff_base_functions = col_data.getFTensor1DiffN<3>(gg, 0);
const int final_bb = diagonal_block ? row_bb + 1 : col_nb_dofs / 3;
int col_bb = 0;
for (; col_bb != final_bb; col_bb++) {
auto t_assemble = get_tensor2(locK, 3 * row_bb, 3 * col_bb);
diffDiff(j, l) =
w * row_diff_base_functions(j) * col_diff_base_functions(l);
t_assemble(i, k) += diffDiff(j, l) * commonData.tD(i, j, k, l);
// Next base function for column
++col_diff_base_functions;
}
++row_diff_base_functions;
}
for (; row_bb != row_nb_base_functions; row_bb++) {
++row_diff_base_functions;
}
}
if (diagonal_block) {
for (int row_bb = 0; row_bb != row_nb_dofs / 3; row_bb++) {
int col_bb = 0;
for (; col_bb != row_bb + 1; col_bb++) {
auto t_assemble = get_tensor2(locK, 3 * row_bb, 3 * col_bb);
auto t_off_side = get_tensor2(locK, 3 * col_bb, 3 * row_bb);
t_off_side(i, k) = t_assemble(k, i);
}
}
}
const int *row_ind = &*row_data.getIndices().begin();
const int *col_ind = &*col_data.getIndices().begin();
Mat B = getFEMethod()->ksp_B != PETSC_NULL ? getFEMethod()->ksp_B
: getFEMethod()->ksp_B;
CHKERR MatSetValues(B, row_nb_dofs, row_ind, col_nb_dofs, col_ind,
if (row_type != col_type || row_side != col_side) {
translocK.resize(col_nb_dofs, row_nb_dofs, false);
noalias(translocK) = trans(locK);
CHKERR MatSetValues(B, col_nb_dofs, col_ind, row_nb_dofs, row_ind,
}
}
};
std::vector<EntityHandle> &mapGaussPts;
std::vector<EntityHandle> &map_gauss_pts,
DataAtIntegrationPts &common_data, BlockData &data)
: VolumeElementForcesAndSourcesCore::UserDataOperator(
commonData(common_data), postProcMesh(post_proc_mesh),
mapGaussPts(map_gauss_pts), dAta(data) {
doVertices = true;
doEdges = false;
doTris = false;
doTets = false;
doPrisms = false;
}
PetscErrorCode doWork(int side, EntityType type,
if (type != MBVERTEX)
PetscFunctionReturn(9);
double def_VAL[9];
bzero(def_VAL, 9 * sizeof(double));
if (dAta.tEts.find(getNumeredEntFiniteElementPtr()->getEnt()) ==
dAta.tEts.end()) {
}
Tag th_stress;
CHKERR postProcMesh.tag_get_handle("STRESS", 9, MB_TYPE_DOUBLE, th_stress,
MB_TAG_CREAT | MB_TAG_SPARSE, def_VAL);
Tag th_strain;
CHKERR postProcMesh.tag_get_handle("STRAIN", 9, MB_TYPE_DOUBLE, th_strain,
MB_TAG_CREAT | MB_TAG_SPARSE, def_VAL);
Tag th_psi;
CHKERR postProcMesh.tag_get_handle("ENERGY", 1, MB_TYPE_DOUBLE, th_psi,
MB_TAG_CREAT | MB_TAG_SPARSE, def_VAL);
const double mu = commonData.mU;
for (int gg = 0; gg != nb_gauss_pts; gg++) {
double psi = 0.5 * p * p + mu * strain(i, j) * strain(i, j);
stress(i, j) = 2 * mu * strain(i, j);
stress(1, 1) -= p;
stress(0, 0) -= p;
stress(2, 2) -= p;
CHKERR postProcMesh.tag_set_data(th_psi, &mapGaussPts[gg], 1, &psi);
CHKERR postProcMesh.tag_set_data(th_strain, &mapGaussPts[gg], 1,
&strain(0, 0));
CHKERR postProcMesh.tag_set_data(th_stress, &mapGaussPts[gg], 1,
&stress(0, 0));
++p;
}
}
};
#endif //__ELASTICITYMIXEDFORMULATION_HPP__
MoFEMFunctionReturnHot
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:447
DataAtIntegrationPts::getParameters
MoFEMErrorCode getParameters()
Definition: ElasticityMixedFormulation.hpp:43
MoFEM::DataOperator::doTris
bool & doTris
\deprectaed
Definition: DataOperators.hpp:93
OpPostProcStress::commonData
DataAtIntegrationPts & commonData
Definition: ElasticityMixedFormulation.hpp:431
MoFEM::MatSetValues
MoFEMErrorCode MatSetValues(Mat M, const EntitiesFieldData::EntData &row_data, const EntitiesFieldData::EntData &col_data, const double *ptr, InsertMode iora)
Assemble PETSc matrix.
Definition: EntitiesFieldData.hpp:1631
OpAssembleK::translocK
MatrixDouble translocK
Definition: ElasticityMixedFormulation.hpp:315
OpAssembleG::OpAssembleG
OpAssembleG(DataAtIntegrationPts &common_data, BlockData &data)
Definition: ElasticityMixedFormulation.hpp:222
FTensor::Tensor1< double, 3 >
EntityHandle
OpAssembleK::dAta
BlockData & dAta
Definition: ElasticityMixedFormulation.hpp:319
BlockData::pOisson
double pOisson
Definition: ElasticityMixedFormulation.hpp:16
DataAtIntegrationPts::mU
double mU
Definition: ElasticityMixedFormulation.hpp:29
BlockData::oRder
int oRder
Definition: gel_analysis.cpp:75
OpPostProcStress::postProcMesh
moab::Interface & postProcMesh
Definition: ElasticityMixedFormulation.hpp:432
OpPostProcStress::OpPostProcStress
OpPostProcStress(moab::Interface &post_proc_mesh, std::vector< EntityHandle > &map_gauss_pts, DataAtIntegrationPts &common_data, BlockData &data)
Definition: ElasticityMixedFormulation.hpp:436
MoFEM::Exceptions::MoFEMErrorCode
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
_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
MoFEM::Types::MatrixDouble
UBlasMatrix< double > MatrixDouble
Definition: Types.hpp:77
MoFEM::DataOperator::doPrisms
bool & doPrisms
\deprectaed
Definition: DataOperators.hpp:95
DataAtIntegrationPts::pOisson
double pOisson
Definition: ElasticityMixedFormulation.hpp:26
DataAtIntegrationPts::setOfBlocksData
std::map< int, BlockData > setOfBlocksData
Definition: ElasticityMixedFormulation.hpp:31
MoFEM::DataOperator::doTets
bool & doTets
\deprectaed
Definition: DataOperators.hpp:94
DataAtIntegrationPts::getBlockData
MoFEMErrorCode getBlockData(BlockData &data)
Definition: ElasticityMixedFormulation.hpp:52
MoFEM::DataOperator::doVertices
bool & doVertices
\deprectaed If false skip vertices
Definition: DataOperators.hpp:90
sdf.r
int r
Definition: sdf.py:8
FTensor::Tensor2< double, 3, 3 >
MoFEM::DeprecatedCoreInterface
Deprecated interface functions.
Definition: DeprecatedCoreInterface.hpp:16
OpAssembleP::doWork
PetscErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
Definition: ElasticityMixedFormulation.hpp:129
MoFEM::Interface
DeprecatedCoreInterface Interface
Definition: Interface.hpp:1975
DataAtIntegrationPts::lAmbda
double lAmbda
Definition: ElasticityMixedFormulation.hpp:28
c
const double c
speed of light (cm/ns)
Definition: initial_diffusion.cpp:39
OpPostProcStress::mapGaussPts
std::vector< EntityHandle > & mapGaussPts
Definition: ElasticityMixedFormulation.hpp:433
OpAssembleP::translocP
MatrixDouble translocP
Definition: ElasticityMixedFormulation.hpp:119
CHKERR
#define CHKERR
Inline error check.
Definition: definitions.h:535
OpAssembleG::locG
MatrixDouble locG
Definition: ElasticityMixedFormulation.hpp:219
OpAssembleK::diffDiff
FTensor::Tensor2< double, 3, 3 > diffDiff
Definition: ElasticityMixedFormulation.hpp:316
OpAssembleG::commonData
DataAtIntegrationPts & commonData
Definition: ElasticityMixedFormulation.hpp:218
MoFEM::CoreInterface::get_moab
virtual moab::Interface & get_moab()=0
OpAssembleK
Assemble K *.
Definition: ElasticityMixedFormulation.hpp:311
OpAssembleP::OpAssembleP
OpAssembleP(DataAtIntegrationPts &common_data, BlockData &data)
Definition: ElasticityMixedFormulation.hpp:122
BlockData
Definition: gel_analysis.cpp:74
convert.type
type
Definition: convert.py:64
OpAssembleG::doWork
PetscErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
Definition: ElasticityMixedFormulation.hpp:229
OpAssembleG::dAta
BlockData & dAta
Definition: ElasticityMixedFormulation.hpp:220
MoFEM::getFTensor0FromVec
static auto getFTensor0FromVec(ublas::vector< T, A > &data)
Get tensor rank 0 (scalar) form data vector.
Definition: Templates.hpp:135
MAT_ELASTICSET
@ MAT_ELASTICSET
block name is "MAT_ELASTIC"
Definition: definitions.h:159
DataAtIntegrationPts::DataAtIntegrationPts
DataAtIntegrationPts()
Definition: HookeElement.hpp:96
DataAtIntegrationPts
Definition: HookeElement.hpp:79
OpPostProcStress::dAta
BlockData & dAta
Definition: ElasticityMixedFormulation.hpp:434
DataAtIntegrationPts::mField
MoFEM::Interface & mField
Definition: ElasticityMixedFormulation.hpp:104
MoFEM::ForcesAndSourcesCore::UserDataOperator::getNumeredEntFiniteElementPtr
boost::shared_ptr< const NumeredEntFiniteElement > getNumeredEntFiniteElementPtr() const
Return raw pointer to NumeredEntFiniteElement.
Definition: ForcesAndSourcesCore.hpp:999
OpAssembleP
Assemble P *.
Definition: ElasticityMixedFormulation.hpp:114
\deprectaed
Definition: DataOperators.hpp:92
BlockData::iD
int iD
Definition: ElasticityMixedFormulation.hpp:13
OpPostProcStress::doWork
PetscErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Definition: ElasticityMixedFormulation.hpp:451
i
FTensor::Index< 'i', SPACE_DIM > i
Definition: hcurl_divergence_operator_2d.cpp:27
EntData
EntitiesFieldData::EntData EntData
Definition: child_and_parent.cpp:37
FTensor::Index< 'i', 3 >
OpAssembleP::commonData
DataAtIntegrationPts & commonData
Definition: ElasticityMixedFormulation.hpp:117
BlockData::BlockData
BlockData()
Definition: gel_analysis.cpp:76
BlockData::tEts
Range tEts
Definition: ElasticityMixedFormulation.hpp:17
DataAtIntegrationPts::yOung
double yOung
Definition: ElasticityMixedFormulation.hpp:27
Range
MoFEM::DataOperator::doEdges
bool & doEdges
\deprectaed If false skip edges
Definition: DataOperators.hpp:91
DataAtIntegrationPts::pPtr
boost::shared_ptr< VectorDouble > pPtr
Definition: ElasticityMixedFormulation.hpp:23
Definition: ElasticityMixedFormulation.hpp:22
OpAssembleK::locK
MatrixDouble locK
Definition: ElasticityMixedFormulation.hpp:314
OpAssembleP::dAta
BlockData & dAta
Definition: ElasticityMixedFormulation.hpp:120
UserDataOperator
ForcesAndSourcesCore::UserDataOperator UserDataOperator
Definition: HookeElement.hpp:75
BLOCKSET
@ BLOCKSET
Definition: definitions.h:148
j
FTensor::Index< 'j', 3 > j
Definition: matrix_function.cpp:19
FTensor::Ddg< double, 3, 3 >
MoFEM::Exceptions::ierr
static MoFEMErrorCodeGeneric< PetscErrorCode > ierr
Definition: Exceptions.hpp:76
mu
double mu
Definition: dynamic_first_order_con_law.cpp:98
lambda
static double lambda
Definition: incompressible_elasticity.cpp:199
DataAtIntegrationPts::tD
FTensor::Ddg< double, 3, 3 > tD
Definition: ElasticityMixedFormulation.hpp:24
DataAtIntegrationPts::setBlocks
MoFEMErrorCode setBlocks()
Definition: ElasticityMixedFormulation.hpp:86
OpAssembleK::commonData
DataAtIntegrationPts & commonData
Definition: ElasticityMixedFormulation.hpp:318
OpAssembleK::doWork
PetscErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
Definition: ElasticityMixedFormulation.hpp:327
MoFEM::Types::VectorDouble
UBlasVector< double > VectorDouble
Definition: Types.hpp:68
sdf_wavy_2d.w
int w
Definition: sdf_wavy_2d.py:7
BlockData::yOung
double yOung
Definition: ElasticityMixedFormulation.hpp:15
m
FTensor::Index< 'm', 3 > m
Definition: shallow_wave.cpp:80
k
FTensor::Index< 'k', 3 > k
Definition: matrix_function.cpp:20
OpAssembleP::locP
MatrixDouble locP
Definition: ElasticityMixedFormulation.hpp:118
MoFEMFunctionReturn
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:416
OpAssembleK::OpAssembleK
OpAssembleK(DataAtIntegrationPts &common_data, BlockData &data, bool symm=true)
Definition: ElasticityMixedFormulation.hpp:321
MoFEMFunctionBegin
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:346
l
FTensor::Index< 'l', 3 > l
Definition: matrix_function.cpp:21
MoFEM::ForcesAndSourcesCore::UserDataOperator::OPROW
@ OPROW
operator doWork function is executed on FE rows
Definition: ForcesAndSourcesCore.hpp:567
OpAssembleG
Assemble G *.
Definition: ElasticityMixedFormulation.hpp:215