13static char help[] =
"...\n\n";
49 GAUSS>::OpGradGradSymTensorGradGrad<1, 1, SPACE_DIM, 0>;
58constexpr double mu = 1;
59constexpr double t = 1;
73 return cos(2 * x * M_PI) * sin(2 * y * M_PI);
82 auto mat_D_ptr = boost::make_shared<MatrixDouble>(
a *
a, 1);
83 auto t_D = getFTensor4DdgFromMat<2, 2, 0>(*(mat_D_ptr));
84 constexpr double t3 =
t *
t *
t;
85 constexpr double A =
mu * t3 / 12;
86 constexpr double B =
lambda * t3 / 12;
96std::array<std::vector<VectorInt>, 2>
98std::array<std::vector<MatrixDouble>, 2>
100std::array<std::vector<MatrixDouble>, 2>
129 boost::shared_ptr<MatrixDouble> d_mat_ptr);
135 boost::shared_ptr<FaceSideEle>
204 auto get_skin = [&]() {
210 MOAB_THROW(skin.find_skin(0, body_ents,
false, skin_ents));
213 skin_ents.merge(verts);
218 auto remove_dofs_from_problem = [&](
Range &&skin) {
222 CHKERR problem_mng->removeDofsOnEntities(
simple->getProblemName(),
"U",
228 CHKERR remove_dofs_from_problem(get_skin());
240 auto rule_lhs = [](int, int,
int p) ->
int {
return 2 *
p; };
241 auto rule_rhs = [](int, int,
int p) ->
int {
return 2 *
p; };
242 auto rule_2 = [
this](int, int, int) {
return 2 *
order; };
244 CHKERR pipeline_mng->setDomainLhsIntegrationRule(rule_lhs);
245 CHKERR pipeline_mng->setDomainRhsIntegrationRule(rule_rhs);
247 CHKERR pipeline_mng->setSkeletonLhsIntegrationRule(rule_2);
248 CHKERR pipeline_mng->setSkeletonRhsIntegrationRule(rule_2);
249 CHKERR pipeline_mng->setBoundaryLhsIntegrationRule(rule_2);
250 CHKERR pipeline_mng->setBoundaryRhsIntegrationRule(rule_2);
252 auto det_ptr = boost::make_shared<VectorDouble>();
253 auto jac_ptr = boost::make_shared<MatrixDouble>();
254 auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
255 auto base_mass_ptr = boost::make_shared<MatrixDouble>();
256 auto data_l2_ptr = boost::make_shared<EntitiesFieldData>(MBENTITYSET);
260 auto push_ho_direcatives = [&](
auto &pipeline) {
264 BaseDerivatives::SecondDerivative, base_mass_ptr, data_l2_ptr,
272 auto push_jacobian = [&](
auto &pipeline) {
283 push_ho_direcatives(pipeline_mng->getOpDomainLhsPipeline());
284 push_jacobian(pipeline_mng->getOpDomainLhsPipeline());
286 pipeline_mng->getOpDomainLhsPipeline().push_back(
291 pipeline_mng->getOpDomainRhsPipeline().push_back(
295 auto side_fe_ptr = boost::make_shared<FaceSideEle>(
mField);
296 push_ho_direcatives(side_fe_ptr->getOpPtrVector());
297 push_jacobian(side_fe_ptr->getOpPtrVector());
301 pipeline_mng->getOpSkeletonLhsPipeline().push_back(
315 auto ksp_solver = pipeline_mng->createKSP();
316 CHKERR KSPSetFromOptions(ksp_solver);
317 CHKERR KSPSetUp(ksp_solver);
320 auto dm =
simple->getDM();
327 CHKERR VecGhostUpdateBegin(
D, INSERT_VALUES, SCATTER_FORWARD);
328 CHKERR VecGhostUpdateEnd(
D, INSERT_VALUES, SCATTER_FORWARD);
341 pipeline_mng->getSkeletonRhsFE().reset();
342 pipeline_mng->getSkeletonLhsFE().reset();
343 pipeline_mng->getBoundaryRhsFE().reset();
344 pipeline_mng->getBoundaryLhsFE().reset();
346 auto post_proc_fe = boost::make_shared<PostProcEle>(
mField);
348 auto u_ptr = boost::make_shared<VectorDouble>();
349 post_proc_fe->getOpPtrVector().push_back(
354 post_proc_fe->getOpPtrVector().push_back(
356 new OpPPMap(post_proc_fe->getPostProcMesh(),
357 post_proc_fe->getMapGaussPts(),
367 pipeline_mng->getDomainRhsFE() = post_proc_fe;
368 CHKERR pipeline_mng->loopFiniteElements();
369 CHKERR post_proc_fe->writeFile(
"out_result.h5m");
391int main(
int argc,
char *argv[]) {
398 auto core_log = logging::core::get();
406 DMType dm_name =
"DMMOFEM";
411 moab::Core mb_instance;
412 moab::Interface &moab = mb_instance;
431 std::string col_field_name)
438 const auto nb_in_loop = getFEMethod()->nInTheLoop;
440 auto clear = [](
auto nb) {
446 if (type == MBVERTEX) {
447 areaMap[nb_in_loop] = getMeasure();
448 senseMap[nb_in_loop] = getEdgeSense();
457 const auto nb_dofs = data.
getIndices().size();
461 data.
getN(BaseDerivatives::FirstDerivative));
463 data.
getN(BaseDerivatives::SecondDerivative));
471 &*base_mat.data().begin());
474template <
typename T>
inline auto get_ntensor(T &base_mat,
int gg,
int bb) {
475 double *ptr = &base_mat(gg, bb);
480 double *ptr = &*base_mat.data().begin();
481 return getFTensor1FromPtr<2>(ptr);
486 double *ptr = &base_mat(gg, 2 * bb);
487 return getFTensor1FromPtr<2>(ptr);
491 double *ptr = &*base_mat.data().begin();
497 double *ptr = &base_mat(gg, 4 * bb);
507 boost::shared_ptr<MatrixDouble> mat_d_ptr)
509 dMatPtr(mat_d_ptr) {}
517 const auto in_the_loop =
526 t_normal.normalize();
530 auto t_D = getFTensor4DdgFromMat<SPACE_DIM, SPACE_DIM, 0>(*
dMatPtr);
533 const size_t nb_integration_pts =
getGaussPts().size2();
537 const double beta =
static_cast<double>(
nitsche) / (in_the_loop + 1);
539 auto integrate = [&](
auto sense_row,
auto &row_ind,
auto &row_diff,
540 auto &row_diff2,
auto sense_col,
auto &col_ind,
541 auto &col_diff,
auto &col_diff2) {
546 const auto nb_rows = row_ind.size();
547 const auto nb_cols = col_ind.size();
549 const auto nb_row_base_functions = row_diff.size2() /
SPACE_DIM;
554 locMat.resize(nb_rows, nb_cols,
false);
564 for (
size_t gg = 0; gg != nb_integration_pts; ++gg) {
569 auto t_mat =
locMat.data().begin();
573 for (; rr != nb_rows; ++rr) {
576 t_mv(
i,
j) = t_D(
i,
j,
k,
l) * t_diff2_row_base(
k,
l);
580 t_vn_plus(
i,
j) = beta * (
phi * t_mv(
i,
j) /
p);
582 t_vn(
i,
j) = (t_diff_row_base(
j) * (t_normal(
i) * sense_row)) -
590 for (
size_t cc = 0; cc != nb_cols; ++cc) {
593 t_mu(
i,
j) = t_D(
i,
j,
k,
l) * t_diff2_col_base(
k,
l);
597 t_un(
i,
j) = -
p * ((t_diff_col_base(
j) * (t_normal(
i) * sense_col) -
598 beta * t_mu(
i,
j) /
p));
601 *t_mat -= alpha * (t_vn(
i,
j) * t_un(
i,
j));
602 *t_mat -= alpha * (t_vn_plus(
i,
j) * (beta * t_mu(
i,
j)));
620 for (; rr < nb_row_base_functions; ++rr) {
629 CHKERR ::MatSetValues(
getKSPB(), nb_rows, &*row_ind.begin(),
630 col_ind.size(), &*col_ind.begin(),
631 &*
locMat.data().begin(), ADD_VALUES);
640 const auto sense_row =
senseMap[s0];
645 const auto sense_col =
senseMap[s1];
void simple(double P1[], double P2[], double P3[], double c[], const int N)
ElementsAndOps< SPACE_DIM >::DomainEle DomainEle
ElementsAndOps< SPACE_DIM >::BoundaryEle BoundaryEle
Kronecker Delta class symmetric.
#define CATCH_ERRORS
Catch errors.
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
#define MOAB_THROW(err)
Check error code of MoAB function and throw MoFEM exception.
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
@ L2
field with C-1 continuity
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define CHKERR
Inline error check.
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
ElementsAndOps< SPACE_DIM >::FaceSideEle FaceSideEle
PetscErrorCode DMoFEMMeshToLocalVector(DM dm, Vec l, InsertMode mode, ScatterMode scatter_mode)
set local (or ghosted) vector values on mesh for partition only
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
auto createDMVector(DM dm)
Get smart vector from DM.
static LoggerType & setLog(const std::string channel)
Set ans resset chanel logger.
#define MOFEM_LOG(channel, severity)
Log.
#define MOFEM_LOG_TAG(channel, tag)
Tag channel.
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpGradGrad< 1, 1, SPACE_DIM > OpDomainGradGrad
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
implementation of Data Operators for Forces and Sources
PetscErrorCode PetscOptionsGetInt(PetscOptions *, const char pre[], const char name[], PetscInt *ivalue, PetscBool *set)
PetscErrorCode PetscOptionsGetScalar(PetscOptions *, const char pre[], const char name[], PetscScalar *dval, PetscBool *set)
SmartPetscObj< Vec > vectorDuplicate(Vec vec)
Create duplicate vector of smart vector.
FTensor::Tensor2_symmetric< FTensor::PackPtr< double *, 4 >, 2 > getFTensor2SymmetricLowerFromPtr< 2 >(double *ptr)
OpPostProcMapInMoab< SPACE_DIM, SPACE_DIM > OpPPMap
constexpr double mu
lame parameter
auto get_diff_ntensor(T &base_mat)
std::array< std::vector< MatrixDouble >, 2 > diff2BaseSideMap
FTensor::Index< 'j', SPACE_DIM > j
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpGradGradSymTensorGradGrad< 1, 1, SPACE_DIM, 0 > OpDomainPlateStiffness
std::array< std::vector< MatrixDouble >, 2 > diffBaseSideMap
FTensor::Index< 'k', SPACE_DIM > k
FTensor::Index< 'i', SPACE_DIM > i
auto get_ntensor(T &base_mat)
constexpr int SPACE_DIM
dimension of space
FTensor::Index< 'l', SPACE_DIM > l
constexpr int FIELD_DIM
dimension of approx. field
auto plate_stiffness
get fourth-order constitutive tensor
std::array< double, 2 > areaMap
std::array< std::vector< VectorInt >, 2 > indicesSideMap
indices on rows for left hand-side
std::array< int, 2 > senseMap
constexpr int BASE_DIM
dimension of base
constexpr double t
plate stiffness
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpSource< BASE_DIM, FIELD_DIM > OpDomainPlateLoad
auto get_diff2_ntensor(T &base_mat)
constexpr auto field_name
virtual moab::Interface & get_moab()=0
static MoFEMErrorCode Initialize(int *argc, char ***args, const char file[], const char help[])
Initializes the MoFEM database PETSc, MOAB and MPI.
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Deprecated interface functions.
Data on single entity (This is passed as argument to DataOperator::doWork)
MatrixDouble & getN(const FieldApproximationBase base)
get base functions this return matrix (nb. of rows is equal to nb. of Gauss pts, nb....
const VectorInt & getIndices() const
Get global indices of dofs on entity.
default operator for TRI element
auto getFTensor1Normal()
get normal as tensor
default operator for Face element
Base face element used to integrate on skeleton.
auto getFTensor0IntegrationWeight()
Get integration weights.
double getMeasure() const
get measure of element
MatrixDouble & getGaussPts()
matrix of integration (Gauss) points for Volume Element
static boost::shared_ptr< SinkType > createSink(boost::shared_ptr< std::ostream > stream_ptr, std::string comm_filter)
Create a sink object.
static boost::shared_ptr< std::ostream > getStrmWorld()
Get the strm world object.
Get value at integration points for scalar field.
Post post-proc data at points from hash maps.
Modify integration weights on face to take in account higher-order geometry.
PipelineManager interface.
boost::shared_ptr< FEMethod > & getDomainLhsFE()
MoFEM::FaceElementForcesAndSourcesCore FaceEle
MoFEM::EdgeElementForcesAndSourcesCore EdgeEle
Problem manager is used to build and partition problems.
Simple interface for fast problem set-up.
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
Operator tp collect data from elements on the side of Edge/Face.
OpCalculateSideData(std::string field_name, std::string col_field_name, UserDataOperator::OpType type)
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
Operator the left hand side matrix.
boost::shared_ptr< MatrixDouble > dMatPtr
boost::shared_ptr< FaceSideEle > sideFEPtr
pointer to element to get data on edge/face sides
MatrixDouble locMat
local operator matrix
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
OpH1LhsSkeleton(boost::shared_ptr< FaceSideEle > side_fe_ptr, boost::shared_ptr< MatrixDouble > d_mat_ptr)
Construct a new OpH1LhsSkeleton.
MoFEMErrorCode setupProblem()
[Read mesh]
MoFEM::Interface & mField
MoFEMErrorCode runProblem()
[Output results]
MoFEMErrorCode boundaryCondition()
[Set up problem]
MoFEMErrorCode outputResults()
[Solve system]
MoFEMErrorCode solveSystem()
[Push operators to pipeline]
MoFEMErrorCode assembleSystem()
[Boundary condition]
Plate(MoFEM::Interface &m_field)
MoFEMErrorCode readMesh()
[Read mesh]