v0.13.0
approx_sphere.cpp
/**
* \file approx_sphere.cpp
* \example approx_sphere.cpp
*
*/
/* This file is part of MoFEM.
* MoFEM is free software: you can redistribute it and/or modify it under
* the terms of the GNU Lesser General Public License as published by the
* Free Software Foundation, either version 3 of the License, or (at your
* option) any later version.
*
* MoFEM is distributed in the hope that it will be useful, but WITHOUT
* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
* FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
* License for more details.
*
* You should have received a copy of the GNU Lesser General Public
* License along with MoFEM. If not, see <http://www.gnu.org/licenses/>. */
#include <MoFEM.hpp>
using namespace MoFEM;
static char help[] = "...\n\n";
#include <BasicFiniteElements.hpp>
constexpr int FM_DIM = 2;
template <int DIM> struct ElementsAndOps {};
template <> struct ElementsAndOps<2> {
};
FormsIntegrators<DomainEleOp>::Assembly<PETSC>::OpBase;
constexpr double a = 1;
constexpr double a2 = a * a;
constexpr double a4 = a2 * a2;
constexpr double A = 6371220;
auto res_J = [](const double x, const double y, const double z) {
const double res = (x * x + y * y + z * z - a2);
return res;
};
auto res_J_dx = [](const double x, const double y, const double z) {
const double res = res_J(x, y, z);
return FTensor::Tensor1<double, 3>{res * (2 * x), res * (2 * y),
res * (2 * z)};
};
auto lhs_J_dx2 = [](const double x, const double y, const double z) {
const double res = res_J(x, y, z);
(res * 2 + (4 * x * x)),
(4 * y * x),
(4 * z * x),
(4 * x * y),
(2 * res + (4 * y * y)),
(4 * z * y),
(4 * x * z),
(4 * y * z),
(2 * res + (4 * z * z))};
};
struct OpRhs : public AssemblyDomainEleOp {
OpRhs(const std::string field_name, boost::shared_ptr<MatrixDouble> x_ptr,
boost::shared_ptr<MatrixDouble> dot_x_ptr)
: AssemblyDomainEleOp(field_name, field_name, AssemblyDomainEleOp::OPROW),
xPtr(x_ptr), xDotPtr(dot_x_ptr) {}
auto t_w = getFTensor0IntegrationWeight();
auto t_row_base = row_data.getFTensor0N();
auto t_x0 = getFTensor1CoordsAtGaussPts();
auto t_x = getFTensor1FromMat<3>(*xPtr);
auto t_dot_x = getFTensor1FromMat<3>(*xDotPtr);
auto t_normal = getFTensor1NormalsAtGaussPts();
for (int gg = 0; gg != nbIntegrationPts; gg++) {
FTensor::Tensor1<double, 3> t_n{t_x0(0), t_x0(1), t_x0(2)};
t_n.normalize();
t_P(i, j) = t_n(i) * t_n(j);
t_Q(i, j) = t_kd(i, j) - t_P(i, j);
auto t_J_res = res_J_dx(t_x(0), t_x(1), t_x(2));
const double alpha = t_w;
auto t_nf = getFTensor1FromArray<3, 3>(locF);
double l = std::sqrt(t_normal(i) * t_normal(i));
t_res(i) =
alpha * l * ((t_P(i, k) * t_J_res(k) + t_Q(i, k) * t_dot_x(k)));
int rr = 0;
for (; rr != nbRows / 3; ++rr) {
t_nf(j) += t_row_base * t_res(j);
++t_row_base;
++t_nf;
}
for (; rr < nbRowBaseFunctions; ++rr) {
++t_row_base;
}
++t_w;
++t_x;
++t_dot_x;
++t_x0;
++t_normal;
}
}
private:
boost::shared_ptr<MatrixDouble> xPtr;
boost::shared_ptr<MatrixDouble> xDotPtr;
};
struct OpLhs : public AssemblyDomainEleOp {
OpLhs(const std::string field_name, boost::shared_ptr<MatrixDouble> x_ptr,
boost::shared_ptr<MatrixDouble> dot_x_ptr)
: AssemblyDomainEleOp(field_name, field_name,
AssemblyDomainEleOp::OPROWCOL),
xPtr(x_ptr), xDotPtr(dot_x_ptr) {
this->sYmm = false;
}
auto t_w = getFTensor0IntegrationWeight();
auto t_row_base = row_data.getFTensor0N();
auto t_x0 = getFTensor1CoordsAtGaussPts();
auto t_x = getFTensor1FromMat<3>(*xPtr);
auto t_dot_x = getFTensor1FromMat<3>(*xDotPtr);
auto t_normal = getFTensor1NormalsAtGaussPts();
auto get_t_mat = [&](const int rr) {
&locMat(rr + 0, 0), &locMat(rr + 0, 1), &locMat(rr + 0, 2),
&locMat(rr + 1, 0), &locMat(rr + 1, 1), &locMat(rr + 1, 2),
&locMat(rr + 2, 0), &locMat(rr + 2, 1), &locMat(rr + 2, 2)};
};
const double ts_a = getTSa();
for (int gg = 0; gg != nbIntegrationPts; gg++) {
FTensor::Tensor1<double, 3> t_n{t_x0(0), t_x0(1), t_x0(2)};
t_n.normalize();
t_P(i, j) = t_n(i) * t_n(j);
t_Q(i, j) = t_kd(i, j) - t_P(i, j);
auto t_J_lhs = lhs_J_dx2(t_x(0), t_x(1), t_x(2));
double l = std::sqrt(t_normal(i) * t_normal(i));
const double alpha = t_w;
t_lhs(i, j) =
(alpha * l) * (t_P(i, k) * t_J_lhs(k, j) + t_Q(i, j) * ts_a);
int rr = 0;
for (; rr != nbRows / 3; rr++) {
auto t_col_base = col_data.getFTensor0N(gg, 0);
auto t_mat = get_t_mat(3 * rr);
for (int cc = 0; cc != nbCols / 3; cc++) {
const double rc = t_row_base * t_col_base;
t_mat(i, j) += rc * t_lhs(i, j);
++t_col_base;
++t_mat;
}
++t_row_base;
}
for (; rr < nbRowBaseFunctions; ++rr)
++t_row_base;
++t_w;
++t_x;
++t_x0;
++t_normal;
}
}
private:
boost::shared_ptr<MatrixDouble> xPtr;
boost::shared_ptr<MatrixDouble> xDotPtr;
};
struct OpError : public DomainEleOp {
OpError(const std::string field_name, boost::shared_ptr<MatrixDouble> x_ptr)
: DomainEleOp(field_name, field_name, AssemblyDomainEleOp::OPROW),
xPtr(x_ptr) {
std::fill(&doEntities[MBEDGE], &doEntities[MBMAXTYPE], false);
}
MoFEMErrorCode doWork(int side, EntityType type, EntData &data) {
auto t_w = getFTensor0IntegrationWeight();
auto t_x = getFTensor1FromMat<3>(*xPtr);
auto t_normal = getFTensor1NormalsAtGaussPts();
auto nb_integration_pts = getGaussPts().size2();
double error = 0;
for (int gg = 0; gg != nb_integration_pts; gg++) {
double l = std::sqrt(t_normal(i) * t_normal(i));
error += t_w * l * std::abs((t_x(i) * t_x(i) - A * A));
++t_w;
++t_x;
++t_normal;
}
CHKERR VecSetValue(errorVec, 0, error, ADD_VALUES);
}
static SmartPetscObj<Vec> errorVec;
private:
boost::shared_ptr<MatrixDouble> xPtr;
};
SmartPetscObj<Vec> OpError::errorVec;
struct ApproxSphere {
ApproxSphere(MoFEM::Interface &m_field) : mField(m_field) {}
MoFEMErrorCode runProblem();
private:
MoFEMErrorCode getOptions();
MoFEMErrorCode readMesh();
MoFEMErrorCode setupProblem();
MoFEMErrorCode setOPs();
MoFEMErrorCode solveSystem();
MoFEMErrorCode outputResults();
};
//! [Run programme]
CHKERR getOptions();
CHKERR readMesh();
CHKERR setupProblem();
CHKERR setOPs();
CHKERR solveSystem();
CHKERR outputResults();
}
//! [Run programme]
}
//! [Read mesh]
auto simple = mField.getInterface<Simple>();
CHKERR simple->getOptions();
CHKERR simple->loadFile();
}
//! [Read mesh]
//! [Set up problem]
auto simple = mField.getInterface<Simple>();
CHKERR simple->addDomainField("HO_POSITIONS", H1, AINSWORTH_LEGENDRE_BASE, 3);
CHKERR simple->addDataField("HO_POSITIONS", H1, AINSWORTH_LEGENDRE_BASE, 3);
int order = 3;
CHKERR PetscOptionsGetInt(PETSC_NULL, "", "-order", &order, PETSC_NULL);
CHKERR simple->setFieldOrder("HO_POSITIONS", order);
CHKERR simple->setUp();
}
//! [Set up problem]
//! [Push operators to pipeline]
PipelineManager *pipeline_mng = mField.getInterface<PipelineManager>();
auto integration_rule = [](int, int, int approx_order) {
return 3 * approx_order;
};
CHKERR pipeline_mng->setDomainRhsIntegrationRule(integration_rule);
CHKERR pipeline_mng->setDomainLhsIntegrationRule(integration_rule);
auto x_ptr = boost::make_shared<MatrixDouble>();
auto dot_x_ptr = boost::make_shared<MatrixDouble>();
auto det_ptr = boost::make_shared<VectorDouble>();
auto jac_ptr = boost::make_shared<MatrixDouble>();
auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
auto def_ops = [&](auto &pipeline) {
pipeline.push_back(
new OpCalculateVectorFieldValues<3>("HO_POSITIONS", x_ptr));
pipeline.push_back(
new OpCalculateVectorFieldValuesDot<3>("HO_POSITIONS", dot_x_ptr));
};
def_ops(pipeline_mng->getOpDomainRhsPipeline());
def_ops(pipeline_mng->getOpDomainLhsPipeline());
pipeline_mng->getOpDomainRhsPipeline().push_back(
new OpRhs("HO_POSITIONS", x_ptr, dot_x_ptr));
pipeline_mng->getOpDomainLhsPipeline().push_back(
new OpLhs("HO_POSITIONS", x_ptr, dot_x_ptr));
}
//! [Push operators to pipeline]
//! [Solve]
// Project HO geometry from mesh
Projection10NodeCoordsOnField ent_method_material(mField, "HO_POSITIONS");
CHKERR mField.loop_dofs("HO_POSITIONS", ent_method_material);
auto *simple = mField.getInterface<Simple>();
auto *pipeline_mng = mField.getInterface<PipelineManager>();
auto dm = simple->getDM();
ts = pipeline_mng->createTSIM();
double ftime = 1;
CHKERR TSSetMaxSteps(ts, 1);
CHKERR TSSetExactFinalTime(ts, TS_EXACTFINALTIME_MATCHSTEP);
auto T = smartCreateDMVector(simple->getDM());
CHKERR DMoFEMMeshToLocalVector(simple->getDM(), T, INSERT_VALUES,
SCATTER_FORWARD);
CHKERR TSSetSolution(ts, T);
CHKERR TSSetFromOptions(ts);
CHKERR TSSolve(ts, NULL);
CHKERR TSGetTime(ts, &ftime);
CHKERR mField.getInterface<FieldBlas>()->fieldScale(A, "HO_POSITIONS");
}
//! [Solve]
auto simple = mField.getInterface<Simple>();
auto dm = simple->getDM();
auto post_proc_fe = boost::make_shared<PostProcFaceOnRefinedMesh>(mField);
post_proc_fe->generateReferenceElementMesh();
post_proc_fe->addFieldValuesPostProc("HO_POSITIONS");
CHKERR DMoFEMLoopFiniteElements(dm, "dFE", post_proc_fe);
CHKERR post_proc_fe->writeFile("out_approx.h5m");
auto error_fe = boost::make_shared<DomainEle>(mField);
auto x_ptr = boost::make_shared<MatrixDouble>();
auto det_ptr = boost::make_shared<VectorDouble>();
auto jac_ptr = boost::make_shared<MatrixDouble>();
auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
error_fe->getOpPtrVector().push_back(
new OpGetHONormalsOnFace("HO_POSITIONS"));
error_fe->getOpPtrVector().push_back(
new OpCalculateVectorFieldValues<3>("HO_POSITIONS", x_ptr));
error_fe->getOpPtrVector().push_back(new OpError("HO_POSITIONS", x_ptr));
error_fe->preProcessHook = [&]() {
MOFEM_LOG("EXAMPLE", Sev::inform) << "Create vec ";
mField.get_comm(), (!mField.get_comm_rank()) ? 1 : 0, 1);
VecZeroEntries(OpError::errorVec);
};
error_fe->postProcessHook = [&]() {
CHKERR VecAssemblyBegin(OpError::errorVec);
CHKERR VecAssemblyEnd(OpError::errorVec);
double error2;
CHKERR VecSum(OpError::errorVec, &error2);
MOFEM_LOG("EXAMPLE", Sev::inform)
<< "Error " << std::sqrt(error2 / (4 * M_PI * A * A));
};
CHKERR DMoFEMLoopFiniteElements(dm, "dFE", error_fe);
CHKERR simple->deleteDM();
CHKERR simple->deleteFiniteElements();
if (mField.get_comm_size() > 1)
CHKERR mField.get_moab().write_file("out_ho_mesh.h5m", "MOAB",
"PARALLEL=WRITE_PART");
else
CHKERR mField.get_moab().write_file("out_ho_mesh.h5m");
}
//! [Postprocess results]
int main(int argc, char *argv[]) {
// Initialisation of MoFEM/PETSc and MOAB data structures
const char param_file[] = "param_file.petsc";
auto core_log = logging::core::get();
core_log->add_sink(
LogManager::setLog("EXAMPLE");
MOFEM_LOG_TAG("EXAMPLE", "example");
try {
//! [Register MoFEM discrete manager in PETSc]
DMType dm_name = "DMMOFEM";
//! [Register MoFEM discrete manager in PETSc
//! [Create MoAB]
moab::Core mb_instance; ///< mesh database
moab::Interface &moab = mb_instance; ///< mesh database interface
//! [Create MoAB]
//! [Create MoFEM]
MoFEM::Core core(moab); ///< finite element database
MoFEM::Interface &m_field = core; ///< finite element database insterface
//! [Create MoFEM]
//! [ApproxSphere]
ApproxSphere ex(m_field);
CHKERR ex.runProblem();
//! [ApproxSphere]
}
}
std::string param_file
void simple(double P1[], double P2[], double P3[], double c[], const int N)
Definition: acoustic.cpp:69
constexpr double a2
int main(int argc, char *argv[])
[Postprocess results]
EntitiesFieldData::EntData EntData
constexpr double a
static char help[]
constexpr double a4
constexpr double A
auto lhs_J_dx2
constexpr int FM_DIM
FTensor::Index< 'j', 3 > j
FTensor::Index< 'i', 3 > i
auto res_J_dx
auto res_J
FTensor::Index< 'k', 3 > k
ElementsAndOps< SPACE_DIM >::DomainEle DomainEle
Kronecker Delta class.
Tensor1< T, Tensor_Dim > normalize()
constexpr double alpha
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:385
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:73
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:460
@ H1
continuous field
Definition: definitions.h:98
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:359
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:429
#define CHKERR
Inline error check.
Definition: definitions.h:548
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:453
auto integration_rule
constexpr auto t_kd
auto smartCreateDMVector
Get smart vector from DM.
Definition: DMMoFEM.hpp:956
PetscErrorCode DMoFEMMeshToLocalVector(DM dm, Vec l, InsertMode mode, ScatterMode scatter_mode)
set local (or ghosted) vector values on mesh for partition only
Definition: DMMMoFEM.cpp:481
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition: DMMMoFEM.cpp:59
PetscErrorCode DMoFEMLoopFiniteElements(DM dm, const char fe_name[], MoFEM::FEMethod *method, CacheTupleWeakPtr cache_ptr=CacheTupleSharedPtr())
Executes FEMethod for finite elements in DM.
Definition: DMMMoFEM.cpp:544
FaceElementForcesAndSourcesCoreSwitch< 0 > FaceElementForcesAndSourcesCore
Face finite element default.
static LoggerType & setLog(const std::string channel)
Set ans resset chanel logger.
Definition: LogManager.cpp:364
#define MOFEM_LOG(channel, severity)
Log.
Definition: LogManager.hpp:311
#define MOFEM_LOG_TAG(channel, tag)
Tag channel.
Definition: LogManager.hpp:342
FTensor::Index< 'l', 3 > l
const double T
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:67
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:21
PetscErrorCode PetscOptionsGetInt(PetscOptions *, const char pre[], const char name[], PetscInt *ivalue, PetscBool *set)
auto createSmartVectorMPI
Create MPI Vector.
CoreTmp< 0 > Core
Definition: Core.hpp:1096
DeprecatedCoreInterface Interface
Definition: Interface.hpp:1965
static constexpr int approx_order
OpBaseImpl< PETSC, EdgeEleOp > OpBase
Definition: radiation.cpp:41
MoFEMErrorCode readMesh()
[Read mesh]
MoFEMErrorCode solveSystem()
[Push operators to pipeline]
MoFEMErrorCode setupProblem()
[Read mesh]
MoFEMErrorCode runProblem()
[Run programme]
MoFEMErrorCode setOPs()
[Set up problem]
MoFEMErrorCode getOptions()
[Run programme]
MoFEMErrorCode outputResults()
[Solve]
Core (interface) class.
Definition: Core.hpp:92
static MoFEMErrorCode Initialize(int *argc, char ***args, const char file[], const char help[])
Initializes the MoFEM database PETSc, MOAB and MPI.
Definition: Core.cpp:85
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition: Core.cpp:125
Deprecated interface functions.
Data on single entity (This is passed as argument to DataOperator::doWork)
static boost::shared_ptr< SinkType > createSink(boost::shared_ptr< std::ostream > stream_ptr, std::string comm_filter)
Create a sink object.
Definition: LogManager.cpp:279
static boost::shared_ptr< std::ostream > getStrmWorld()
Get the strm world object.
Definition: LogManager.cpp:323
VolumeElementForcesAndSourcesCoreBase::UserDataOperator UserDataOperator
static SmartPetscObj< Vec > errorVec
Postprocess on face.