9#ifndef __ELASTICITYMIXEDFORMULATION_HPP__
10#define __ELASTICITYMIXEDFORMULATION_HPP__
23 boost::shared_ptr<VectorDouble>
pPtr;
40 CHKERRABORT(PETSC_COMM_WORLD,
ierr);
45 PetscOptionsBegin(PETSC_COMM_WORLD,
"",
"Problem",
"none");
66 tD(0, 0, 0, 0) = 2 *
mU;
73 tD(1, 1, 1, 1) = 2 *
mU;
80 tD(2, 2, 2, 2) = 2 *
mU;
90 CHKERR it->getAttributeDataStructure(mydata);
91 int id = it->getMeshsetId();
114 :
public VolumeElementForcesAndSourcesCore::UserDataOperator {
128 PetscErrorCode
doWork(
int row_side,
int col_side, EntityType row_type,
130 EntitiesFieldData::EntData &row_data,
131 EntitiesFieldData::EntData &col_data) {
134 const int row_nb_dofs = row_data.getIndices().size();
137 const int col_nb_dofs = col_data.getIndices().size();
141 if (
dAta.
tEts.find(getNumeredEntFiniteElementPtr()->getEnt()) ==
147 locP.resize(row_nb_dofs, col_nb_dofs,
false);
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();
164 if (coefficient != 0.) {
165 for (
int gg = 0; gg != row_nb_gauss_pts; gg++) {
168 double w = getVolume() * getGaussPts()(3, gg);
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++) {
176 locP(row_bb, col_bb) -=
177 w * row_base_functions * col_base_functions * coefficient;
179 ++col_base_functions;
181 ++row_base_functions;
183 for (; row_bb != row_nb_base_functions; row_bb++) {
184 ++row_base_functions;
191 getFEMethod()->ksp_B, row_nb_dofs, &*row_data.getIndices().begin(),
192 col_nb_dofs, &*col_data.getIndices().begin(), &
locP(0, 0), ADD_VALUES);
195 if (row_side != col_side || row_type != col_type) {
196 translocP.resize(col_nb_dofs, row_nb_dofs,
false);
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),
215 :
public VolumeElementForcesAndSourcesCore::UserDataOperator {
228 PetscErrorCode
doWork(
int row_side,
int col_side, EntityType row_type,
230 EntitiesFieldData::EntData &row_data,
231 EntitiesFieldData::EntData &col_data) {
235 const int row_nb_dofs = row_data.getIndices().size();
238 const int col_nb_dofs = col_data.getIndices().size();
242 if (
dAta.
tEts.find(getNumeredEntFiniteElementPtr()->getEnt()) ==
249 locG.resize(row_nb_dofs, col_nb_dofs,
false);
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>();
262 for (
int gg = 0; gg != row_nb_gauss_pts; gg++) {
265 double w = getVolume() * getGaussPts()(3, gg);
268 for (; row_bb != row_nb_dofs / 3; row_bb++) {
270 t1(
i) = w * row_diff_base_functions(
i);
272 auto base_functions = col_data.getFTensor0N(gg, 0);
273 for (
int col_bb = 0; col_bb != col_nb_dofs; col_bb++) {
276 &
locG(3 * row_bb + 1, col_bb),
277 &
locG(3 * row_bb + 2, col_bb));
279 k(
i) += t1(
i) * base_functions;
282 ++row_diff_base_functions;
284 for (; row_bb != row_nb_base_functions; row_bb++) {
285 ++row_diff_base_functions;
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(),
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(),
311 :
public VolumeElementForcesAndSourcesCore::UserDataOperator {
326 PetscErrorCode
doWork(
int row_side,
int col_side, EntityType row_type,
328 EntitiesFieldData::EntData &row_data,
329 EntitiesFieldData::EntData &col_data) {
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));
339 const int row_nb_dofs = row_data.getIndices().size();
342 const int col_nb_dofs = col_data.getIndices().size();
345 if (
dAta.
tEts.find(getNumeredEntFiniteElementPtr()->getEnt()) ==
351 const bool diagonal_block =
352 (row_type == col_type) && (row_side == col_side);
355 locK.resize(row_nb_dofs, col_nb_dofs,
false);
358 const int row_nb_gauss_pts = row_data.getN().size1();
359 const int row_nb_base_functions = row_data.getN().size2();
361 auto row_diff_base_functions = row_data.getFTensor1DiffN<3>();
369 for (
int gg = 0; gg != row_nb_gauss_pts; gg++) {
372 double w = getVolume() * getGaussPts()(3, gg);
375 for (; row_bb != row_nb_dofs / 3; row_bb++) {
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;
380 for (; col_bb != final_bb; col_bb++) {
382 auto t_assemble = get_tensor2(
locK, 3 * row_bb, 3 * col_bb);
385 w * row_diff_base_functions(
j) * col_diff_base_functions(
l);
389 ++col_diff_base_functions;
392 ++row_diff_base_functions;
394 for (; row_bb != row_nb_base_functions; row_bb++) {
395 ++row_diff_base_functions;
399 if (diagonal_block) {
400 for (
int row_bb = 0; row_bb != row_nb_dofs / 3; row_bb++) {
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);
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);
417 if (row_type != col_type || row_side != col_side) {
418 translocK.resize(col_nb_dofs, row_nb_dofs,
false);
420 CHKERR MatSetValues(
B, col_nb_dofs, col_ind, row_nb_dofs, row_ind,
436 std::vector<EntityHandle> &map_gauss_pts,
450 PetscErrorCode
doWork(
int side, EntityType type,
451 EntitiesFieldData::EntData &data) {
453 if (type != MBVERTEX)
454 PetscFunctionReturn(9);
456 bzero(def_VAL, 9 *
sizeof(
double));
466 MB_TAG_CREAT | MB_TAG_SPARSE, def_VAL);
469 MB_TAG_CREAT | MB_TAG_SPARSE, def_VAL);
472 MB_TAG_CREAT | MB_TAG_SPARSE, def_VAL);
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);
489 stress(
i,
j) = 2 *
mu * strain(
i,
j);
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"
#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
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
UBlasVector< double > VectorDouble
FTensor::Index< 'm', 3 > m
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
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
moab::Interface & postProcMesh
PetscErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)