static char help[] =
"...\n\n";
#include <BasicFiniteElements.hpp>
};
};
double operator()(const double x, const double y, const double z) {
return sin(x * 10.) * cos(y * 10.);
}
};
using OpDomainMass = FormsIntegrators<DomainEleOp>::Assembly<
private:
Simple *simpleInterface;
boost::shared_ptr<VectorDouble> approxVals;
SmartPetscObj<Vec> L2Vec;
SmartPetscObj<Vec> resVec;
};
boost::shared_ptr<CommonData> commonDataPtr;
template <
int FIELD_DIM>
struct OpError;
};
boost::shared_ptr<CommonData> commonDataPtr;
OpError(boost::shared_ptr<CommonData> &common_data_ptr)
if (
const size_t nb_dofs = data.
getIndices().size()) {
const int nb_integration_pts = getGaussPts().size2();
auto t_w = getFTensor0IntegrationWeight();
auto t_coords = getFTensor1CoordsAtGaussPts();
nf.clear();
const double volume = getMeasure();
double error = 0;
for (int gg = 0; gg != nb_integration_pts; ++gg) {
const double alpha = t_w * volume;
t_coords(2));
error +=
alpha * pow(diff, 2);
for (
size_t r = 0;
r != nb_dofs; ++
r) {
nf[
r] +=
alpha * t_row_base * diff;
++t_row_base;
}
++t_w;
++t_val;
++t_coords;
}
CHKERR VecSetValue(commonDataPtr->L2Vec,
index, error, ADD_VALUES);
}
}
};
}
CHKERR mField.getInterface(simpleInterface);
CHKERR simpleInterface->getOptions();
CHKERR simpleInterface->loadFile();
}
CHKERR simpleInterface->setUp();
}
auto rule = [](
int,
int,
int p) ->
int {
return 2 *
p; };
PipelineManager *pipeline_mng = mField.getInterface<PipelineManager>();
CHKERR pipeline_mng->setDomainLhsIntegrationRule(rule);
CHKERR pipeline_mng->setDomainRhsIntegrationRule(rule);
}
commonDataPtr = boost::make_shared<CommonData>();
mField.get_comm(), (!mField.get_comm_rank()) ? 1 : 0, 1);
commonDataPtr->approxVals = boost::make_shared<VectorDouble>();
}
PipelineManager *pipeline_mng = mField.getInterface<PipelineManager>();
auto beta = [](const double, const double, const double) { return 1; };
pipeline_mng->getOpDomainLhsPipeline().push_back(
pipeline_mng->getOpDomainRhsPipeline().push_back(
}
PipelineManager *pipeline_mng = mField.getInterface<PipelineManager>();
auto solver = pipeline_mng->createKSP();
CHKERR KSPSetFromOptions(solver);
auto dm = simpleInterface->getDM();
CHKERR VecGhostUpdateBegin(
D, INSERT_VALUES, SCATTER_FORWARD);
CHKERR VecGhostUpdateEnd(
D, INSERT_VALUES, SCATTER_FORWARD);
}
PipelineManager *pipeline_mng = mField.getInterface<PipelineManager>();
pipeline_mng->getDomainLhsFE().reset();
auto post_proc_fe = boost::make_shared<PostProcFaceOnRefinedMesh>(mField);
post_proc_fe->generateReferenceElementMesh();
pipeline_mng->getDomainRhsFE() = post_proc_fe;
CHKERR pipeline_mng->loopFiniteElements();
CHKERR post_proc_fe->writeFile(
"out_approx.h5m");
}
PipelineManager *pipeline_mng = mField.getInterface<PipelineManager>();
pipeline_mng->getDomainLhsFE().reset();
pipeline_mng->getDomainRhsFE().reset();
pipeline_mng->getOpDomainRhsPipeline().clear();
pipeline_mng->getOpDomainRhsPipeline().push_back(
new OpCalculateScalarFieldValues(
FIELD_NAME, commonDataPtr->approxVals));
pipeline_mng->getOpDomainRhsPipeline().push_back(
CHKERR pipeline_mng->loopFiniteElements();
CHKERR VecAssemblyBegin(commonDataPtr->L2Vec);
CHKERR VecAssemblyEnd(commonDataPtr->L2Vec);
CHKERR VecAssemblyBegin(commonDataPtr->resVec);
CHKERR VecAssemblyEnd(commonDataPtr->resVec);
double nrm2;
CHKERR VecNorm(commonDataPtr->resVec, NORM_2, &nrm2);
const double *array;
CHKERR VecGetArrayRead(commonDataPtr->L2Vec, &array);
if (mField.get_comm_rank() == 0)
PetscPrintf(PETSC_COMM_SELF, "Error %6.4e Vec norm %6.4e\n", std::sqrt(array[0]),
nrm2);
CHKERR VecRestoreArrayRead(commonDataPtr->L2Vec, &array);
constexpr
double eps = 1e-8;
"Not converged solution");
}
int main(
int argc,
char *argv[]) {
try {
DMType dm_name = "DMMOFEM";
}
}
int main(int argc, char *argv[])
EntitiesFieldData::EntData EntData
constexpr char FIELD_NAME[]
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.
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
MoFEMErrorCode boundaryCondition()
[Set up problem]
MoFEMErrorCode assembleSystem()
[Boundary condition]
MoFEMErrorCode readMesh()
[Run problem]
MoFEMErrorCode setIntegrationRules()
[Set up problem]
static ApproxFieldFunction< FIELD_DIM > approxFunction
MoFEMErrorCode checkResults()
[Postprocess results]
MoFEMErrorCode solveSystem()
[Solve]
MoFEMErrorCode createCommonData()
[Create common data]
MoFEMErrorCode runProblem()
[Create common data]
MoFEMErrorCode setupProblem()
[Read mesh]
MoFEMErrorCode outputResults()
[Solve]
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.
MoFEM::FaceElementForcesAndSourcesCore FaceEle
VolumeElementForcesAndSourcesCoreBase::UserDataOperator UserDataOperator