static char help[] =
"...\n\n";
FaceElementForcesAndSourcesCore::NO_HO_GEOMETRY |
FaceElementForcesAndSourcesCore::NO_CONTRAVARIANT_TRANSFORM_HDIV |
FaceElementForcesAndSourcesCore::NO_COVARIANT_TRANSFORM_HCURL>;
static double fUn(const double x, const double y) {
for (
int i = 0;
i <=
o; ++
i) {
r += pow(x,
i) * pow(y,
j);
}
}
}
for (
int i = 0;
i <=
o; ++
i) {
r(0) +=
i > 0 ?
i * pow(x,
i - 1) * pow(y,
j) : 0;
r(1) +=
j > 0 ?
j * pow(x,
i) * pow(y,
j - 1) : 0;
}
}
}
}
};
sYmm = false;
}
MoFEMErrorCode doWork(
int row_side,
int col_side, EntityType row_type,
EntityType col_type,
EntData &row_data,
const int nb_dofs_row = row_data.
getIndices().size();
if (nb_dofs_row == 0)
const int nb_dofs_col = col_data.
getIndices().size();
if (nb_dofs_col == 0)
nA.resize(nb_dofs_row, nb_dofs_col, false);
nA.clear();
const int nb_gauss_pts = row_data.
getN().size1();
for (int gg = 0; gg != nb_gauss_pts; gg++) {
const double val = getArea() * getGaussPts()(2, gg);
for (int rr = 0; rr != nb_dofs_row; rr++) {
for (int cc = 0; cc != nb_dofs_col; cc++) {
nA(rr, cc) += val * t_base_row * t_base_col;
++t_base_col;
}
++t_base_row;
}
}
&*nA.data().begin(), ADD_VALUES);
}
};
const int nb_dofs = data.getIndices().size();
if (nb_dofs == 0)
const int nb_gauss_pts = data.getN().size1();
nF.resize(nb_dofs, false);
nF.clear();
auto t_base_fun = data.getFTensor0N();
for (int gg = 0; gg != nb_gauss_pts; gg++) {
const double val = getArea() * getGaussPts()(2, gg);
for (int bb = 0; bb != nb_dofs; bb++) {
const double x = getCoordsAtGaussPts()(gg, 0);
const double y = getCoordsAtGaussPts()(gg, 1);
++t_base_fun;
}
}
ADD_VALUES);
}
};
const bool checkGradients;
:
FaceEleOp(
"FIELD1", OPROW), vAls(vals), diffVals(diff_vals),
checkGradients(check_grads) {}
const int nb_gauss_pts = getGaussPts().size2();
vAls.resize(nb_gauss_pts, false);
diffVals.resize(2, nb_gauss_pts, false);
vAls.clear();
diffVals.clear();
}
const int nb_dofs = data.getIndices().size();
if (nb_dofs) {
auto t_base_fun = data.getFTensor0N();
for (int gg = 0; gg != nb_gauss_pts; gg++) {
auto t_data = data.getFTensor0FieldData();
for (int bb = 0; bb != nb_dofs; bb++) {
t_vals += t_base_fun * t_data;
++t_base_fun;
++t_data;
}
++t_vals;
}
if (checkGradients) {
auto t_diff_vals = getFTensor1FromMat<2>(diffVals);
auto t_diff_base_fun = data.getFTensor1DiffN<2>();
for (int gg = 0; gg != nb_gauss_pts; gg++) {
auto t_data = data.getFTensor0FieldData();
for (int bb = 0; bb != nb_dofs; bb++) {
t_diff_vals(
i) += t_diff_base_fun(
i) * t_data;
++t_diff_base_fun;
++t_data;
}
++t_diff_vals;
}
}
}
}
};
boost::shared_ptr<VectorDouble> ptrVals;
boost::shared_ptr<MatrixDouble> ptrDiffVals;
const bool checkGradients;
boost::shared_ptr<VectorDouble> &ptr_vals,
boost::shared_ptr<MatrixDouble> &ptr_diff_vals,
bool check_grads)
:
FaceEleOp(
"FIELD1", OPROW), vAls(vals), diffVals(diff_vals),
ptrVals(ptr_vals), ptrDiffVals(ptr_diff_vals),
checkGradients(check_grads) {
std::fill(&doEntities[MBEDGE], &doEntities[MBMAXTYPE], false);
}
const int nb_gauss_pts = data.getN().size1();
for (int gg = 0; gg != nb_gauss_pts; gg++) {
const double x = getCoordsAtGaussPts()(gg, 0);
const double y = getCoordsAtGaussPts()(gg, 1);
double err_val = std::fabs(delta_val * delta_val);
cerr << err_val << " : " << t_vals << endl;
err_val);
err_val = std::abs(t_vals - t_ptr_vals);
"Wrong value from operator %4.3e", err_val);
++t_vals;
++t_ptr_vals;
}
if (checkGradients) {
auto t_diff_vals = getFTensor1FromMat<2>(diffVals);
auto t_ptr_diff_vals = getFTensor1FromMat<2>(*ptrDiffVals);
for (int gg = 0; gg != nb_gauss_pts; gg++) {
const double x = getCoordsAtGaussPts()(gg, 0);
const double y = getCoordsAtGaussPts()(gg, 1);
t_delta_diff_val(
i) = t_diff_vals(
i) - t_diff_anal(
i);
double err_diff_val = sqrt(t_delta_diff_val(
i) * t_delta_diff_val(
i));
cerr << err_diff_val <<
" : " << sqrt(t_diff_vals(
i) * t_diff_vals(
i))
<< " : " << t_diff_vals(0) / t_diff_anal(0) << " "
<< t_diff_vals(1) / t_diff_anal(1) << endl;
"Wrong derivative of value %4.3e", err_diff_val);
t_delta_diff_val(
i) = t_diff_vals(
i) - t_ptr_diff_vals(
i);
err_diff_val = sqrt(t_delta_diff_val(
i) * t_delta_diff_val(
i));
"Wrong direvatives from operator %4.3e", err_diff_val);
++t_diff_vals;
++t_ptr_diff_vals;
}
}
}
};
int main(
int argc,
char *argv[]) {
try {
DMType dm_name = "DMMOFEM";
PipelineManager *pipeline_mng = m_field.
getInterface<PipelineManager>();
CHKERR simple_interface->getOptions();
CHKERR simple_interface->loadFile(
"",
"");
enum bases {
AINSWORTH,
AINSWORTH_LOBATTO,
DEMKOWICZ,
BERNSTEIN,
LASBASETOP
};
const char *list_bases[] = {"ainsworth", "ainsworth_labatto", "demkowicz",
"bernstein"};
PetscBool flg;
PetscInt choice_base_value = AINSWORTH;
LASBASETOP, &choice_base_value, &flg);
if (flg != PETSC_TRUE)
if (choice_base_value == AINSWORTH)
if (choice_base_value == AINSWORTH_LOBATTO)
else if (choice_base_value == DEMKOWICZ)
else if (choice_base_value == BERNSTEIN)
enum spaces { H1SPACE, L2SPACE, LASBASETSPACE };
const char *list_spaces[] = {"h1", "l2"};
PetscInt choice_space_value = H1SPACE;
LASBASETSPACE, &choice_space_value, &flg);
if (flg != PETSC_TRUE)
if (choice_space_value == H1SPACE)
else if (choice_space_value == L2SPACE)
CHKERR simple_interface->addDomainField(
"FIELD1", space, base, 1);
CHKERR simple_interface->setUp();
auto dm = simple_interface->getDM();
auto assemble_matrices_and_vectors = [&]() {
pipeline_mng->getOpDomainRhsPipeline().push_back(
new OpMakeHighOrderGeometryWeightsOnFace());
pipeline_mng->getOpDomainRhsPipeline().push_back(
new OpAssembleVec());
pipeline_mng->getOpDomainLhsPipeline().push_back(
pipeline_mng->getOpDomainRhsPipeline().push_back(
new OpMakeHighOrderGeometryWeightsOnFace());
pipeline_mng->getOpDomainLhsPipeline().push_back(
new OpAssembleMat());
auto integration_rule = [](
int,
int,
int p_data) {
return 2 * p_data + 1;
};
CHKERR pipeline_mng->setDomainRhsIntegrationRule(integration_rule);
CHKERR pipeline_mng->setDomainLhsIntegrationRule(integration_rule);
};
auto solve_problem = [&] {
auto solver = pipeline_mng->createKSP();
CHKERR KSPSetFromOptions(solver);
CHKERR VecGhostUpdateBegin(
D, INSERT_VALUES, SCATTER_FORWARD);
CHKERR VecGhostUpdateEnd(
D, INSERT_VALUES, SCATTER_FORWARD);
};
auto check_solution = [&] {
boost::shared_ptr<VectorDouble> ptr_values(
new VectorDouble());
boost::shared_ptr<MatrixDouble> ptr_diff_vals(
new MatrixDouble());
pipeline_mng->getDomainLhsFE().reset();
pipeline_mng->getOpDomainRhsPipeline().clear();
pipeline_mng->getOpDomainRhsPipeline().push_back(
pipeline_mng->getOpDomainRhsPipeline().push_back(
new OpSetInvJacH1ForFace(inv_jac));
pipeline_mng->getOpDomainRhsPipeline().push_back(
pipeline_mng->getOpDomainRhsPipeline().push_back(
new OpCalculateScalarFieldValues("FIELD1", ptr_values));
if (choice_space_value == H1SPACE) {
pipeline_mng->getOpDomainRhsPipeline().push_back(
new OpCalculateScalarFieldGradient<2>("FIELD1", ptr_diff_vals));
}
pipeline_mng->getOpDomainRhsPipeline().push_back(
choice_space_value == H1SPACE));
CHKERR pipeline_mng->loopFiniteElements();
};
CHKERR assemble_matrices_and_vectors();
}
}
ForcesAndSourcesCore::UserDataOperator UserDataOperator
MoFEM::FaceElementForcesAndSourcesCore FaceEle
#define CATCH_ERRORS
Catch errors.
FieldApproximationBase
approximation base
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
@ AINSWORTH_BERNSTEIN_BEZIER_BASE
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
FieldSpace
approximation spaces
@ L2
field with C-1 continuity
#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.
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.
FTensor::Index< 'j', 3 > j
FTensor::Index< 'i', 3 > i
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
ublas::vector< double, DoubleAllocator > VectorDouble
ublas::matrix< double, ublas::row_major, DoubleAllocator > MatrixDouble
implementation of Data Operators for Forces and Sources
MoFEMErrorCode MatSetValues(Mat M, const DataForcesAndSourcesCore::EntData &row_data, const DataForcesAndSourcesCore::EntData &col_data, const double *ptr, InsertMode iora)
Assemble PETSc matrix.
OpCalculateInvJacForFaceImpl< 2 > OpCalculateInvJacForFace
static FTensor::Tensor0< FTensor::PackPtr< double *, 1 > > 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)
DeprecatedCoreInterface Interface
SmartPetscObj< Vec > smartVectorDuplicate(SmartPetscObj< Vec > &vec)
Create duplicate vector of smart vector.
MoFEMErrorCode VecSetValues(Vec V, const DataForcesAndSourcesCore::EntData &data, const double *ptr, InsertMode iora)
Assemble PETSc vector.
const double D
diffusivity
const double r
rate factor
int main(int argc, char *argv[])
static constexpr int approx_order
DataForcesAndSourcesCore::EntData EntData
static FTensor::Tensor1< double, 3 > fUn(const double x, const double y)
static FTensor::Tensor2< double, 3, 2 > diffFun(const double x, const double y)
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.
Data on single entity (This is passed as argument to DataOperator::doWork)
const VectorInt & getIndices() const
Get global indices of dofs on entity.
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 > > getFTensor0N(const FieldApproximationBase base)
Get base function as Tensor0.
Deprecated interface functions.
default operator for TRI element
Face finite element switched.
MoFEMErrorCode getInterface(const MOFEMuuid &uuid, IFACE *&iface) const
Get interface by uuid and return reference to pointer of interface.