26 using namespace MoFEM;
28 static char help[] =
"...\n\n";
30 #include <BasicFiniteElements.hpp>
57 double operator()(
const double x,
const double y,
const double z) {
58 return sin(x * 10.) * cos(y * 10.);
94 boost::shared_ptr<CommonData> commonDataPtr;
104 OpError(boost::shared_ptr<CommonData> &common_data_ptr)
109 if (
const size_t nb_dofs = data.
getIndices().size()) {
111 const int nb_integration_pts = getGaussPts().size2();
112 auto t_w = getFTensor0IntegrationWeight();
114 auto t_coords = getFTensor1CoordsAtGaussPts();
120 const double volume = getMeasure();
124 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
126 const double alpha = t_w * volume;
129 error +=
alpha * pow(diff, 2);
131 for (
size_t r = 0;
r != nb_dofs; ++
r) {
132 nf[
r] +=
alpha * t_row_base * diff;
142 CHKERR VecSetValue(commonDataPtr->L2Vec,
index, error, ADD_VALUES);
155 CHKERR setIntegrationRules();
156 CHKERR createCommonData();
157 CHKERR boundaryCondition();
170 CHKERR mField.getInterface(simpleInterface);
171 CHKERR simpleInterface->getOptions();
172 CHKERR simpleInterface->loadFile();
184 constexpr
int order = 4;
186 CHKERR simpleInterface->setUp();
195 auto rule = [](
int,
int,
int p) ->
int {
return 2 *
p; };
208 commonDataPtr = boost::make_shared<CommonData>();
211 mField.get_comm(), (!mField.get_comm_rank()) ? 1 : 0, 1);
212 commonDataPtr->approxVals = boost::make_shared<VectorDouble>();
225 auto beta = [](
const double,
const double,
const double) {
return 1; };
239 CHKERR KSPSetFromOptions(solver);
242 auto dm = simpleInterface->getDM();
247 CHKERR VecGhostUpdateBegin(
D, INSERT_VALUES, SCATTER_FORWARD);
248 CHKERR VecGhostUpdateEnd(
D, INSERT_VALUES, SCATTER_FORWARD);
258 auto post_proc_fe = boost::make_shared<PostProcFaceOnRefinedMesh>(mField);
259 post_proc_fe->generateReferenceElementMesh();
260 post_proc_fe->addFieldValuesPostProc(
FIELD_NAME);
263 CHKERR post_proc_fe->writeFile(
"out_approx.h5m");
280 CHKERR VecAssemblyBegin(commonDataPtr->L2Vec);
281 CHKERR VecAssemblyEnd(commonDataPtr->L2Vec);
282 CHKERR VecAssemblyBegin(commonDataPtr->resVec);
283 CHKERR VecAssemblyEnd(commonDataPtr->resVec);
285 CHKERR VecNorm(commonDataPtr->resVec, NORM_2, &nrm2);
287 CHKERR VecGetArrayRead(commonDataPtr->L2Vec, &array);
288 if (mField.get_comm_rank() == 0)
289 PetscPrintf(PETSC_COMM_SELF,
"Error %6.4e Vec norm %6.4e\n", std::sqrt(array[0]),
291 CHKERR VecRestoreArrayRead(commonDataPtr->L2Vec, &array);
292 constexpr
double eps = 1e-8;
295 "Not converged solution");
300 int main(
int argc,
char *argv[]) {
309 DMType dm_name =
"DMMOFEM";
int main(int argc, char *argv[])
[Check results]
EntitiesFieldData::EntData EntData
constexpr char FIELD_NAME[]
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpSource< 1, FIELD_DIM > OpDomainSource
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, FIELD_DIM > OpDomainMass
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpSource< 1, FIELD_DIM > OpDomainSource
ElementsAndOps< SPACE_DIM >::DomainEle DomainEle
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, FIELD_DIM > OpDomainMass
#define CATCH_ERRORS
Catch errors.
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
@ MOFEM_DATA_INCONSISTENCY
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define CHKERR
Inline error check.
auto smartCreateDMVector
Get smart vector from DM.
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.
MoFEMErrorCode loopFiniteElements(SmartPetscObj< DM > dm=nullptr)
Iterate finite elements.
boost::ptr_vector< UserDataOperator > & getOpDomainLhsPipeline()
Get the Op Domain Lhs Pipeline object.
SmartPetscObj< KSP > createKSP(SmartPetscObj< DM > dm=nullptr)
Create KSP (linear) solver.
boost::ptr_vector< UserDataOperator > & getOpDomainRhsPipeline()
Get the Op Domain Rhs Pipeline object.
VolumeElementForcesAndSourcesCoreSwitch< 0 > VolumeElementForcesAndSourcesCore
Volume finite element default.
FTensor::Index< 'i', SPACE_DIM > i
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
UBlasVector< double > VectorDouble
implementation of Data Operators for Forces and Sources
auto createSmartVectorMPI
Create MPI Vector.
static auto getFTensor0FromVec(ublas::vector< T, A > &data)
Get tensor rank 0 (scalar) form data vector.
DeprecatedCoreInterface Interface
SmartPetscObj< Vec > smartVectorDuplicate(SmartPetscObj< Vec > &vec)
Create duplicate vector of smart vector.
MoFEMErrorCode VecSetValues(Vec V, const EntitiesFieldData::EntData &data, const double *ptr, InsertMode iora)
Assemble PETSc vector.
const double D
diffusivity
const double r
rate factor
double operator()(const double x, const double y, const double z)
PipelineManager::FaceEle DomainEle
SmartPetscObj< Vec > L2Vec
SmartPetscObj< Vec > resVec
boost::shared_ptr< VectorDouble > approxVals
OpError(boost::shared_ptr< CommonData > &common_data_ptr)
boost::shared_ptr< CommonData > commonDataPtr
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
Operator for linear form, usually to calculate values on right hand side.
MoFEMErrorCode boundaryCondition()
MoFEMErrorCode assembleSystem()
MoFEMErrorCode readMesh()
MoFEMErrorCode setIntegrationRules()
static ApproxFieldFunction< FIELD_DIM > approxFunction
MoFEMErrorCode checkResults()
MoFEMErrorCode solveSystem()
MoFEMErrorCode createCommonData()
Example(MoFEM::Interface &m_field)
MoFEMErrorCode runProblem()
MoFEMErrorCode setupProblem()
MoFEMErrorCode outputResults()
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)
FTensor::Tensor0< FTensor::PackPtr< double *, 1 > > getFTensor0N(const FieldApproximationBase base)
Get base function as Tensor0.
const VectorInt & getIndices() const
Get global indices of dofs on entity.
Get value at integration points for scalar field.
PipelineManager interface.
boost::shared_ptr< FEMethod > & getDomainRhsFE()
boost::shared_ptr< FEMethod > & getDomainLhsFE()
MoFEM::FaceElementForcesAndSourcesCore FaceEle
MoFEMErrorCode setDomainRhsIntegrationRule(RuleHookFun rule)
MoFEMErrorCode setDomainLhsIntegrationRule(RuleHookFun rule)
Simple interface for fast problem set-up.
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
Volume finite element with switches.
VolumeElementForcesAndSourcesCoreBase::UserDataOperator UserDataOperator