Implementation Kirchhoff-Love plate using Discointinous Galerkin (DG-Nitsche method)
static char help[] =
"...\n\n";
};
GAUSS>::OpGradGradSymTensorGradGrad<1, 1, SPACE_DIM, 0>;
-1;
auto source = [](
const double x,
const double y,
const double) {
return cos(2 * x * M_PI) * sin(2 * y * M_PI);
};
auto mat_D_ptr = boost::make_shared<MatrixDouble>(a * a, 1);
auto t_D = getFTensor4DdgFromMat<2, 2, 0>(*(mat_D_ptr));
constexpr double t3 =
t *
t *
t;
constexpr double A =
mu * t3 / 12;
constexpr double B =
lambda * t3 / 12;
return mat_D_ptr;
};
std::array<std::vector<VectorInt>, 2>
std::array<std::vector<MatrixDouble>, 2>
std::array<std::vector<MatrixDouble>, 2>
};
public:
boost::shared_ptr<MatrixDouble> d_mat_ptr);
private:
boost::shared_ptr<FaceSideEle>
boost::shared_ptr<MatrixDouble>
dMatPtr;
};
private:
};
}
PETSC_NULL);
PETSC_NULL);
}
auto get_skin = [&]() {
MOAB_THROW(skin.find_skin(0, body_ents,
false, skin_ents));
skin_ents.merge(verts);
return skin_ents;
};
auto remove_dofs_from_problem = [&](
Range &&skin) {
CHKERR problem_mng->removeDofsOnEntities(
simple->getProblemName(),
"U",
skin, 0, 1);
};
CHKERR remove_dofs_from_problem(get_skin());
}
auto rule_lhs = [](
int,
int,
int p) ->
int {
return 2 *
p; };
auto rule_rhs = [](
int,
int,
int p) ->
int {
return 2 *
p; };
CHKERR pipeline_mng->setDomainLhsIntegrationRule(rule_lhs);
CHKERR pipeline_mng->setDomainRhsIntegrationRule(rule_rhs);
CHKERR pipeline_mng->setSkeletonLhsIntegrationRule(rule_2);
CHKERR pipeline_mng->setSkeletonRhsIntegrationRule(rule_2);
CHKERR pipeline_mng->setBoundaryLhsIntegrationRule(rule_2);
CHKERR pipeline_mng->setBoundaryRhsIntegrationRule(rule_2);
auto det_ptr = boost::make_shared<VectorDouble>();
auto jac_ptr = boost::make_shared<MatrixDouble>();
auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
auto base_mass_ptr = boost::make_shared<MatrixDouble>();
auto data_l2_ptr = boost::make_shared<EntitiesFieldData>(MBENTITYSET);
auto push_ho_direcatives = [&](auto &pipeline) {
BaseDerivatives::SecondDerivative, base_mass_ptr, data_l2_ptr,
};
auto push_jacobian = [&](auto &pipeline) {
pipeline.push_back(
};
push_ho_direcatives(pipeline_mng->getOpDomainLhsPipeline());
push_jacobian(pipeline_mng->getOpDomainLhsPipeline());
pipeline_mng->getOpDomainLhsPipeline().push_back(
pipeline_mng->getOpDomainRhsPipeline().push_back(
auto side_fe_ptr = boost::make_shared<FaceSideEle>(
mField);
push_ho_direcatives(side_fe_ptr->getOpPtrVector());
push_jacobian(side_fe_ptr->getOpPtrVector());
pipeline_mng->getOpSkeletonLhsPipeline().push_back(
}
auto ksp_solver = pipeline_mng->createKSP();
CHKERR KSPSetFromOptions(ksp_solver);
CHKERR VecGhostUpdateBegin(
D, INSERT_VALUES, SCATTER_FORWARD);
CHKERR VecGhostUpdateEnd(
D, INSERT_VALUES, SCATTER_FORWARD);
}
pipeline_mng->getSkeletonRhsFE().reset();
pipeline_mng->getSkeletonLhsFE().reset();
pipeline_mng->getBoundaryRhsFE().reset();
pipeline_mng->getBoundaryLhsFE().reset();
auto post_proc_fe = boost::make_shared<PostProcEle>(
mField);
auto u_ptr = boost::make_shared<VectorDouble>();
post_proc_fe->getOpPtrVector().push_back(
post_proc_fe->getOpPtrVector().push_back(
new OpPPMap(post_proc_fe->getPostProcMesh(),
post_proc_fe->getMapGaussPts(),
{{"U", u_ptr}},
{}, {}, {}
)
);
pipeline_mng->getDomainRhsFE() = post_proc_fe;
CHKERR pipeline_mng->loopFiniteElements();
CHKERR post_proc_fe->writeFile(
"out_result.h5m");
}
}
int main(
int argc,
char *argv[]) {
auto core_log = logging::core::get();
core_log->add_sink(LogManager::createSink(LogManager::getStrmWorld(), "PL"));
LogManager::setLog("PL");
try {
DMType dm_name = "DMMOFEM";
moab::Core mb_instance;
moab::Interface &moab = mb_instance;
}
}
std::string col_field_name)
const auto nb_in_loop = getFEMethod()->nInTheLoop;
auto clear = [](auto nb) {
};
if (type == MBVERTEX) {
areaMap[nb_in_loop] = getMeasure();
if (!nb_in_loop) {
clear(0);
clear(1);
}
}
if (nb_dofs) {
data.
getN(BaseDerivatives::FirstDerivative));
data.
getN(BaseDerivatives::SecondDerivative));
}
}
template <
typename T>
inline auto get_ntensor(
T &base_mat) {
&*base_mat.data().begin());
};
template <
typename T>
inline auto get_ntensor(
T &base_mat,
int gg,
int bb) {
double *ptr = &base_mat(gg, bb);
};
double *ptr = &*base_mat.data().begin();
return getFTensor1FromPtr<2>(ptr);
};
template <typename T>
double *ptr = &base_mat(gg, 2 * bb);
return getFTensor1FromPtr<2>(ptr);
};
double *ptr = &*base_mat.data().begin();
return getFTensor2SymmetricLowerFromPtr<2>(ptr);
};
template <typename T>
double *ptr = &base_mat(gg, 4 * bb);
return getFTensor2SymmetricLowerFromPtr<2>(ptr);
};
boost::shared_ptr<MatrixDouble> mat_d_ptr)
dMatPtr(mat_d_ptr) {}
const auto in_the_loop =
t_normal.normalize();
auto t_D = getFTensor4DdgFromMat<SPACE_DIM, SPACE_DIM, 0>(*
dMatPtr);
const size_t nb_integration_pts =
getGaussPts().size2();
const double beta =
static_cast<double>(
nitsche) / (in_the_loop + 1);
auto integrate = [&](auto sense_row, auto &row_ind, auto &row_diff,
auto &row_diff2, auto sense_col, auto &col_ind,
auto &col_diff, auto &col_diff2) {
const auto nb_rows = row_ind.size();
const auto nb_cols = col_ind.size();
const auto nb_row_base_functions = row_diff.size2() /
SPACE_DIM;
if (nb_cols) {
locMat.resize(nb_rows, nb_cols,
false);
for (size_t gg = 0; gg != nb_integration_pts; ++gg) {
auto t_mat =
locMat.data().begin();
size_t rr = 0;
for (; rr != nb_rows; ++rr) {
t_mv(
i,
j) = t_D(
i,
j,
k,
l) * t_diff2_row_base(
k,
l);
t_vn_plus(
i,
j) = beta * (
phi * t_mv(
i,
j) /
p);
t_vn(
i,
j) = (t_diff_row_base(
j) * (t_normal(
i) * sense_row)) -
for (size_t cc = 0; cc != nb_cols; ++cc) {
t_mu(
i,
j) = t_D(
i,
j,
k,
l) * t_diff2_col_base(
k,
l);
t_un(
i,
j) = -
p * ((t_diff_col_base(
j) * (t_normal(
i) * sense_col) -
*t_mat -=
alpha * (t_vn_plus(
i,
j) * (beta * t_mu(
i,
j)));
++t_diff_col_base;
++t_diff2_col_base;
++t_mat;
}
++t_diff_row_base;
++t_diff2_row_base;
}
for (; rr < nb_row_base_functions; ++rr) {
++t_diff_row_base;
++t_diff2_row_base;
}
++t_w;
}
CHKERR ::MatSetValues(
getKSPB(), nb_rows, &*row_ind.begin(),
col_ind.size(), &*col_ind.begin(),
&*
locMat.data().begin(), ADD_VALUES);
}
};
);
}
}
}
}
}
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 ...
auto get_diff_ntensor(T &base_mat)
auto get_ntensor(T &base_mat)
std::array< double, 2 > areaMap
std::array< int, 2 > senseMap
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.
#define MOFEM_LOG(channel, severity)
Log.
#define MOFEM_LOG_TAG(channel, tag)
Tag channel.
FTensor::Index< 'i', SPACE_DIM > i
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpGradGrad< 1, 1, SPACE_DIM > OpDomainGradGrad
FTensor::Index< 'l', 3 > l
FTensor::Index< 'j', 3 > j
FTensor::Index< 'k', 3 > k
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
UBlasMatrix< double > MatrixDouble
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.
OpPostProcMapInMoab< SPACE_DIM, SPACE_DIM > OpPPMap
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
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 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
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
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()
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]
MoFEMErrorCode readMesh()
[Read mesh]