13static char help[] =
"...\n\n";
33 static double fUn(
const double x,
const double y,
double z) {
36 for (
int i = 0;
i <=
o; ++
i) {
39 r += pow(x,
i) * pow(y,
j);
49 for (
int i = 0;
i <=
o; ++
i) {
52 r(0) +=
i > 0 ?
i * pow(x,
i - 1) * pow(y,
j) : 0;
53 r(1) +=
j > 0 ?
j * pow(x,
i) * pow(y,
j - 1) : 0;
62 static double fUn(
const double x,
const double y,
double z) {
65 for (
int i = 0;
i <=
o; ++
i) {
66 for (
int j = 0;
j <=
o -
i;
j++) {
69 r += pow(x,
i) * pow(y,
j) * pow(z,
k);
81 for (
int i = 0;
i <=
o; ++
i) {
82 for (
int j = 0;
j <=
o -
i;
j++) {
85 r(0) +=
i > 0 ?
i * pow(x,
i - 1) * pow(y,
j) * pow(z,
k) : 0;
86 r(1) +=
j > 0 ?
j * pow(x,
i) * pow(y,
j - 1) * pow(z,
k) : 0;
87 r(2) +=
k > 0 ?
k * pow(x,
i) * pow(y,
j) * pow(z,
k - 1) : 0;
112 if (type == MBVERTEX) {
113 vAls.resize(nb_gauss_pts,
false);
123 <<
"Type " << moab::CN::EntityTypeName(type) <<
" side " << side;
129 for (
int gg = 0; gg != nb_gauss_pts; gg++) {
131 for (
int bb = 0; bb != nb_dofs; bb++) {
132 t_vals += t_base_fun * t_data;
136 const double v = t_vals;
137 if (!std::isnormal(
v))
143 auto t_diff_vals = getFTensor1FromMat<SPACE_DIM>(
diffVals);
145 for (
int gg = 0; gg != nb_gauss_pts; gg++) {
147 for (
int bb = 0; bb != nb_dofs; bb++) {
148 t_diff_vals(
i) += t_diff_base_fun(
i) * t_data;
153 if (!std::isnormal(t_diff_vals(d)))
171 boost::shared_ptr<VectorDouble> &ptr_vals,
172 boost::shared_ptr<MatrixDouble> &ptr_diff_vals,
185 const double eps = 1e-6;
186 const int nb_gauss_pts = data.
getN().size1();
192 for (
int gg = 0; gg != nb_gauss_pts; gg++) {
197 err_val = std::abs(t_vals - t_ptr_vals);
198 MOFEM_LOG(
"AT", Sev::noisy) <<
"Val op error " << err_val;
202 "Wrong value from operator %4.3e", err_val);
209 const double delta_val = t_vals - ApproxFunctions::fUn(x, y, z);
211 err_val = std::fabs(delta_val * delta_val);
212 MOFEM_LOG(
"AT", Sev::verbose) << err_val <<
" : " << t_vals;
223 auto t_diff_vals = getFTensor1FromMat<SPACE_DIM>(
diffVals);
224 auto t_ptr_diff_vals = getFTensor1FromMat<SPACE_DIM>(*
ptrDiffVals);
226 for (
int gg = 0; gg != nb_gauss_pts; gg++) {
231 t_delta_diff_val(
i) = t_diff_vals(
i) - t_ptr_diff_vals(
i);
232 err_diff_val = sqrt(t_delta_diff_val(
i) * t_delta_diff_val(
i));
233 MOFEM_LOG(
"AT", Sev::noisy) <<
"Diff val op error " << err_diff_val;
235 if (err_diff_val >
eps)
237 "Wrong derivatives from operator %4.3e", err_diff_val);
244 auto t_diff_anal = ApproxFunctions::diffFun(x, y, z);
245 t_delta_diff_val(
i) = t_diff_vals(
i) - t_diff_anal(
i);
247 err_diff_val = sqrt(t_delta_diff_val(
i) * t_delta_diff_val(
i));
250 <<
"Diff val " << err_diff_val <<
" : "
251 << sqrt(t_diff_vals(
i) * t_diff_vals(
i)) <<
" : "
252 << t_diff_vals(0) <<
" (" << t_diff_anal(0) <<
") "
253 << t_diff_vals(1) <<
" (" << t_diff_anal(1) <<
") "
254 << t_diff_vals(2) <<
" (" << t_diff_anal(2) <<
")";
257 <<
"Diff val " << err_diff_val <<
" : "
258 << sqrt(t_diff_vals(
i) * t_diff_vals(
i)) <<
" : "
259 << t_diff_vals(0) <<
" (" << t_diff_anal(0) <<
") "
260 << t_diff_vals(1) <<
" (" << t_diff_anal(1) <<
")";
269 MOFEM_LOG(
"AT", Sev::verbose) <<
"Diff val error " << err_diff_val;
270 if (err_diff_val >
eps)
272 "Wrong derivative of value %4.3e %4.3e", err_diff_val,
284int main(
int argc,
char *argv[]) {
290 DMType dm_name =
"DMMOFEM";
293 moab::Core mb_instance;
294 moab::Interface &moab = mb_instance;
297 auto core_log = logging::core::get();
311 simple->getAddBoundaryFE() =
true;
323 const char *list_bases[] = {
"ainsworth",
"ainsworth_labatto",
"demkowicz",
326 PetscInt choice_base_value = AINSWORTH;
328 LASBASETOP, &choice_base_value, &flg);
330 if (flg != PETSC_TRUE)
333 if (choice_base_value == AINSWORTH)
335 if (choice_base_value == AINSWORTH_LOBATTO)
337 else if (choice_base_value == DEMKOWICZ)
339 else if (choice_base_value == BERNSTEIN)
342 enum spaces { H1SPACE, L2SPACE, LASBASETSPACE };
343 const char *list_spaces[] = {
"h1",
"l2"};
344 PetscInt choice_space_value = H1SPACE;
346 LASBASETSPACE, &choice_space_value, &flg);
347 if (flg != PETSC_TRUE)
350 if (choice_space_value == H1SPACE)
352 else if (choice_space_value == L2SPACE)
358 CHKERR simple->addDomainField(
"FIELD1", space, base, 1);
362 CHKERR moab.get_entities_by_dimension(0, 1, edges);
363 CHKERR moab.get_entities_by_dimension(0, 2, faces);
365 if (choice_base_value != BERNSTEIN) {
369 for (
auto e : edges) {
371 rise_order.insert(e);
376 for (
auto f : faces) {
378 rise_order.insert(f);
384 for (
auto e : rise_order) {
386 rise_twice.insert(e);
398 auto volume_adj = [](moab::Interface &moab,
const Field &field,
400 std::vector<EntityHandle> &adjacency) {
403 switch (field.getSpace()) {
405 CHKERR moab.get_connectivity(&fe_ent, 1, adjacency,
true);
407 CHKERR moab.get_adjacencies(&fe_ent, 1, 1,
false, adjacency,
408 moab::Interface::UNION);
411 CHKERR moab.get_adjacencies(&fe_ent, 1, 2,
false, adjacency,
412 moab::Interface::UNION);
413 adjacency.push_back(fe_ent);
415 for (
auto ent : adjacency)
416 fe.getSideNumberPtr(ent);
419 CHKERR moab.get_entities_by_handle(field.getMeshset(), adjacency,
421 for (
auto ent : adjacency) {
424 boost::shared_ptr<SideNumber>(
new SideNumber(ent, -1, 0, 0)));
429 "this field is not implemented for TRI finite element");
436 simple->getDomainFEName(), MBTET, volume_adj);
438 simple->getDomainFEName(), MBHEX, volume_adj);
445 auto dm =
simple->getDM();
448 auto jac_ptr = boost::make_shared<MatrixDouble>();
449 auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
450 auto det_ptr = boost::make_shared<VectorDouble>();
453 auto assemble_matrices_and_vectors = [&]() {
460 pipeline_mng->getOpDomainLhsPipeline(), {NOSPACE});
462 pipeline_mng->getOpDomainRhsPipeline(), {NOSPACE});
467 pipeline_mng->getOpDomainLhsPipeline().push_back(
469 pipeline_mng->getOpDomainLhsPipeline().push_back(
new OpMass(
470 "FIELD1",
"FIELD1", [](
double,
double,
double) {
return 1.; }));
471 pipeline_mng->getOpDomainLhsPipeline().push_back(
474 pipeline_mng->getOpDomainRhsPipeline().push_back(
475 new OpSource(
"FIELD1", ApproxFunctions::fUn));
478 return 2 * p_data + 1;
486 auto solve_problem = [&] {
488 auto solver = pipeline_mng->createKSP();
489 CHKERR KSPSetFromOptions(solver);
492 auto dm =
simple->getDM();
497 CHKERR VecGhostUpdateBegin(
D, INSERT_VALUES, SCATTER_FORWARD);
498 CHKERR VecGhostUpdateEnd(
D, INSERT_VALUES, SCATTER_FORWARD);
504 auto check_solution = [&] {
507 auto ptr_values = boost::make_shared<VectorDouble>();
508 auto ptr_diff_vals = boost::make_shared<MatrixDouble>();
510 pipeline_mng->getDomainLhsFE().reset();
511 pipeline_mng->getOpDomainRhsPipeline().clear();
514 pipeline_mng->getOpDomainRhsPipeline(), {space});
516 pipeline_mng->getOpDomainRhsPipeline().push_back(
518 pipeline_mng->getOpDomainRhsPipeline().push_back(
520 pipeline_mng->getOpDomainRhsPipeline().push_back(
524 vals, diff_vals, ptr_values, ptr_diff_vals,
true));
526 CHKERR pipeline_mng->loopFiniteElements();
531 auto post_proc = [&] {
536 auto get_domain_post_proc_fe = [&](
auto post_proc_mesh_ptr) {
538 boost::make_shared<PostProcBrokenMeshInMoabBaseCont<DomainEle>>(
539 m_field, post_proc_mesh_ptr);
542 post_proc_fe->getOpPtrVector(), {space});
544 auto ptr_values = boost::make_shared<VectorDouble>();
545 auto ptr_diff_vals = boost::make_shared<MatrixDouble>();
547 post_proc_fe->getOpPtrVector().push_back(
549 post_proc_fe->getOpPtrVector().push_back(
553 post_proc_fe->getOpPtrVector().push_back(
557 post_proc_fe->getPostProcMesh(), post_proc_fe->getMapGaussPts(),
559 {{
"FIELD1", ptr_values}},
561 {{
"FIELD1_GRAD", ptr_diff_vals}},
569 auto get_bdy_post_proc_fe = [&](
auto post_proc_mesh_ptr) {
570 auto bdy_post_proc_fe =
571 boost::make_shared<PostProcBrokenMeshInMoabBaseCont<BoundaryEle>>(
572 m_field, post_proc_mesh_ptr);
577 ForcesAndSourcesCore::UserDataOperator::AdjCache>());
579 auto ptr_values = boost::make_shared<VectorDouble>();
580 auto ptr_diff_vals = boost::make_shared<MatrixDouble>();
584 op_loop_side->getOpPtrVector(), {space});
585 op_loop_side->getOpPtrVector().push_back(
587 op_loop_side->getOpPtrVector().push_back(
591 bdy_post_proc_fe->getOpPtrVector().push_back(op_loop_side);
593 bdy_post_proc_fe->getOpPtrVector().push_back(
597 bdy_post_proc_fe->getPostProcMesh(),
598 bdy_post_proc_fe->getMapGaussPts(),
600 {{
"FIELD1", ptr_values}},
602 {{
"FIELD1_GRAD", ptr_diff_vals}},
608 return bdy_post_proc_fe;
611 auto post_proc_mesh_ptr = boost::make_shared<moab::Core>();
612 auto post_proc_begin =
613 boost::make_shared<PostProcBrokenMeshInMoabBaseBegin>(
614 m_field, post_proc_mesh_ptr);
616 boost::make_shared<PostProcBrokenMeshInMoabBaseEnd>(
617 m_field, post_proc_mesh_ptr);
618 auto domain_post_proc_fe = get_domain_post_proc_fe(post_proc_mesh_ptr);
619 auto bdy_post_proc_fe = get_bdy_post_proc_fe(post_proc_mesh_ptr);
623 domain_post_proc_fe);
627 CHKERR post_proc_end->writeFile(
"out_post_proc.h5m");
632 CHKERR assemble_matrices_and_vectors();
void simple(double P1[], double P2[], double P3[], double c[], const int N)
ElementsAndOps< SPACE_DIM >::DomainEle DomainEle
ElementsAndOps< SPACE_DIM >::BoundaryEle BoundaryEle
#define CATCH_ERRORS
Catch errors.
FieldApproximationBase
approximation base
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
@ AINSWORTH_BERNSTEIN_BEZIER_BASE
FieldSpace
approximation spaces
@ L2
field with C-1 continuity
@ NOFIELD
scalar or vector of scalars describe (no true field)
@ HCURL
field with continuous tangents
@ HDIV
field with continuous normal traction
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
@ MOFEM_ATOM_TEST_INVALID
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define CHKERR
Inline error check.
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, SPACE_DIM > OpMass
PetscErrorCode DMoFEMPostProcessFiniteElements(DM dm, MoFEM::FEMethod *method)
execute finite element method for each element in dm (problem)
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.
PetscErrorCode DMoFEMLoopFiniteElements(DM dm, const char fe_name[], MoFEM::FEMethod *method, CacheTupleWeakPtr cache_ptr=CacheTupleSharedPtr())
Executes FEMethod for finite elements in DM.
auto createDMVector(DM dm)
Get smart vector from DM.
PetscErrorCode DMoFEMPreProcessFiniteElements(DM dm, MoFEM::FEMethod *method)
execute finite element method for each element in dm (problem)
multi_index_container< boost::shared_ptr< SideNumber >, indexed_by< ordered_unique< member< SideNumber, EntityHandle, &SideNumber::ent > >, ordered_non_unique< composite_key< SideNumber, const_mem_fun< SideNumber, EntityType, &SideNumber::getEntType >, member< SideNumber, signed char, &SideNumber::side_number > > > > > SideNumber_multiIndex
SideNumber_multiIndex for SideNumber.
virtual MoFEMErrorCode modify_finite_element_adjacency_table(const std::string &fe_name, const EntityType type, ElementAdjacencyFunct function)=0
modify finite element table, only for advanced user
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.
FTensor::Index< 'i', SPACE_DIM > i
const double v
phase velocity of light in medium (cm/ns)
FTensor::Index< 'j', 3 > j
FTensor::Index< 'k', 3 > k
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)
SmartPetscObj< Vec > vectorDuplicate(Vec vec)
Create duplicate vector of smart vector.
static auto getFTensor0FromVec(ublas::vector< T, A > &data)
Get tensor rank 0 (scalar) form data vector.
PetscErrorCode PetscOptionsGetEList(PetscOptions *, const char pre[], const char name[], const char *const *list, PetscInt next, PetscInt *value, PetscBool *set)
OpPostProcMapInMoab< SPACE_DIM, SPACE_DIM > OpPPMap
#define EXECUTABLE_DIMENSION
PipelineManager::ElementsAndOpsByDim< SPACE_DIM >::FaceSideEle EleOnSide
static FTensor::Tensor1< double, 2 > diffFun(const double x, const double y, double z)
static double fUn(const double x, const double y, double z)
static FTensor::Tensor1< double, 3 > diffFun(const double x, const double y, double z)
static double fUn(const double x, const double y, double z)
Add operators pushing bases from local to physical configuration.
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.
std::array< bool, MBMAXTYPE > doEntities
If true operator is executed for entity.
Deprecated interface functions.
Finite element data for entity.
Data on single entity (This is passed as argument to DataOperator::doWork)
FTensor::Tensor1< FTensor::PackPtr< double *, Tensor_Dim >, Tensor_Dim > getFTensor1DiffN(const FieldApproximationBase base)
Get derivatives of base functions.
FTensor::Tensor0< FTensor::PackPtr< double *, 1 > > getFTensor0N(const FieldApproximationBase base)
Get base function as Tensor0.
MatrixDouble & getDiffN(const FieldApproximationBase base)
get derivatives of base functions
MatrixDouble & getN(const FieldApproximationBase base)
get base functions this return matrix (nb. of rows is equal to nb. of Gauss pts, nb....
FTensor::Tensor0< FTensor::PackPtr< double *, 1 > > getFTensor0FieldData()
Resturn scalar files as a FTensor of rank 0.
const VectorInt & getIndices() const
Get global indices of dofs on entity.
VectorDouble & getCoords()
get triangle coordinates
Provide data structure for (tensor) field approximation.
MatrixDouble & getCoordsAtGaussPts()
Gauss points and weight, matrix (nb. of points x 3)
@ OPROW
operator doWork function is executed on FE rows
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 field gradients at integration pts for scalar filed rank 0, i.e. vector field.
Get value at integration points for scalar field.
Element used to execute operators on side of the element.
Post post-proc data at points from hash maps.
Clear Schur complement internal data.
Assemble Schur complement.
PipelineManager interface.
keeps information about side number for the finite element
Simple interface for fast problem set-up.
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
FTensor::Index< 'i', 3 > i
OpCheckValsDiffVals(VectorDouble &vals, MatrixDouble &diff_vals, boost::shared_ptr< VectorDouble > &ptr_vals, boost::shared_ptr< MatrixDouble > &ptr_diff_vals, bool check_grads)
boost::shared_ptr< VectorDouble > ptrVals
const bool checkGradients
FTensor::Index< 'i', SPACE_DIM > i
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
boost::shared_ptr< MatrixDouble > ptrDiffVals
boost::shared_ptr< MatrixDouble > ptrVals
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
const bool checkGradients
OpValsDiffVals(VectorDouble &vals, MatrixDouble &diff_vals, bool check_grads)
FTensor::Index< 'i', SPACE_DIM > i