static char help[] =
"...\n\n";
FaceElementForcesAndSourcesCore::NO_HO_GEOMETRY |
FaceElementForcesAndSourcesCore::NO_CONTRAVARIANT_TRANSFORM_HDIV |
FaceElementForcesAndSourcesCore::NO_COVARIANT_TRANSFORM_HCURL>;
constexpr
double a0 = 0.0;
constexpr
double a1 = 2.0;
constexpr
double a2 = -15.0 *
a0;
constexpr
double a3 = -20.0 / 6 *
a1;
constexpr
double a4 = 15.0 *
a0;
constexpr
double a5 =
a1;
constexpr
double a6 = -
a0;
6 *
a6 * std::pow(x, 5) * std::pow(y, 0) +
5 *
a5 * std::pow(x, 4) * std::pow(y, 1) +
4 *
a4 * std::pow(x, 3) * std::pow(y, 2) +
3 *
a3 * std::pow(x, 2) * std::pow(y, 3) +
2 *
a2 * std::pow(x, 1) * std::pow(y, 4) +
1 *
a1 * std::pow(x, 0) * std::pow(y, 5),
1 *
a5 * std::pow(x, 5) * std::pow(y, 0) +
2 *
a4 * std::pow(x, 4) * std::pow(y, 1) +
3 *
a3 * std::pow(x, 3) * std::pow(y, 2) +
4 *
a2 * std::pow(x, 2) * std::pow(y, 3) +
5 *
a1 * std::pow(x, 1) * std::pow(y, 4) +
6 *
a0 * std::pow(x, 0) * std::pow(y, 5),
0.);
}
const double y) {
30 *
a6 * pow(x, 4) * pow(y, 0) + 20 *
a5 * pow(x, 3) * pow(y, 1) +
12 *
a4 * pow(x, 2) * pow(y, 2) + 6 *
a3 * pow(x, 1) * pow(y, 3) +
2 *
a2 * pow(x, 0) * pow(y, 4),
5 *
a5 * pow(x, 4) * pow(y, 0) + 8 *
a4 * pow(x, 3) * pow(y, 1) +
9 *
a3 * pow(x, 2) * pow(y, 2) + 8 *
a2 * pow(x, 1) * pow(y, 3) +
5 *
a1 * pow(x, 0) * pow(y, 4),
5 *
a5 * pow(x, 4) * pow(y, 0) + 8 *
a4 * pow(x, 3) * pow(y, 1) +
9 *
a3 * pow(x, 2) * pow(y, 2) + 8 *
a2 * pow(x, 1) * pow(y, 3) +
5 *
a1 * pow(x, 0) * pow(y, 4),
2 *
a4 * pow(x, 4) * pow(y, 0) + 6 *
a3 * pow(x, 3) * pow(y, 1) +
12 *
a2 * pow(x, 2) * pow(y, 2) + 20 *
a1 * pow(x, 1) * pow(y, 3) +
30 *
a0 * pow(x, 0) * pow(y, 4),
0., 0.);
}
};
sYmm = false;
}
MoFEMErrorCode doWork(
int row_side,
int col_side, EntityType row_type,
EntityType col_type,
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();
auto t_base_row = row_data.getFTensor1N<3>();
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++) {
auto t_base_col = col_data.getFTensor1N<3>(gg, 0);
for (int cc = 0; cc != nb_dofs_col; cc++) {
nA(rr, cc) += val * t_base_row(
i) * t_base_col(
i);
++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.getFTensor1N<3>();
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);
}
};
:
FaceEleOp(
"FIELD1", OPROW), vAls(vals), diffVals(diff_vals) {}
const int nb_dofs = data.getIndices().size();
if (nb_dofs == 0)
const int nb_gauss_pts = data.getN().size1();
if (
type == MBEDGE && side == 0) {
vAls.resize(3, nb_gauss_pts, false);
diffVals.resize(6, nb_gauss_pts, false);
vAls.clear();
diffVals.clear();
}
auto t_vals = getFTensor1FromMat<3>(vAls);
auto t_diff_vals = getFTensor2FromMat<3, 2>(diffVals);
auto t_base_fun = data.getFTensor1N<3>();
auto t_diff_base_fun = data.getFTensor2DiffN<3, 2>();
for (int gg = 0; gg != nb_gauss_pts; gg++) {
auto t_data = data.getFTensor0FieldData();
for (int bb = 0; bb != nb_dofs; bb++) {
t_vals(
i) += t_base_fun(
i) * t_data;
t_diff_vals(
i,
j) += t_diff_base_fun(
i,
j) * t_data;
++t_base_fun;
++t_diff_base_fun;
++t_data;
}
++t_vals;
++t_diff_vals;
}
}
};
boost::shared_ptr<MatrixDouble> ptrVals;
boost::shared_ptr<VectorDouble> ptrDiv;
boost::shared_ptr<MatrixDouble> ptr_vals,
boost::shared_ptr<VectorDouble> ptr_div)
:
FaceEleOp(
"FIELD1", OPROW), vAls(vals), diffVals(diff_vals),
ptrVals(ptr_vals), ptrDiv(ptr_div) {}
if (
type == MBEDGE && side == 0) {
const int nb_gauss_pts = data.getN().size1();
auto t_vals = getFTensor1FromMat<3>(vAls);
auto t_diff_vals = getFTensor2FromMat<3, 2>(diffVals);
auto t_vals_from_op = getFTensor1FromMat<3>(*ptrVals);
for (int gg = 0; gg != nb_gauss_pts; gg++) {
const double x = getCoordsAtGaussPts()(gg, 0);
const double y = getCoordsAtGaussPts()(gg, 1);
double err_val = sqrt(delta_val(
i) * delta_val(
i));
double err_diff_val = sqrt(delta_diff_val(
i,
j) * delta_diff_val(
i,
j));
"Wrong value %4.3e", err_val);
"Wrong derivative of value %4.3e", err_diff_val);
delta_val(
i) = t_vals(
i) - t_vals_from_op(
i);
err_val = sqrt(delta_val(
i) * delta_val(
i));
"Wrong value from operator %4.3e", err_val);
double div = t_diff_vals(0, 0) + t_diff_vals(1, 1);
double err_div = div - t_div_from_op;
"Wrong divergence from operator %4.3e (%4.3e != %4.3e)",
err_div, div, t_div_from_op);
++t_vals;
++t_diff_vals;
++t_vals_from_op;
++t_div_from_op;
}
}
}
};
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(
"",
"rectangle_tri.h5m");
enum bases { AINSWORTH, DEMKOWICZ, LASBASETOP };
const char *list_bases[] = {"ainsworth", "demkowicz"};
PetscBool flg;
PetscInt choice_base_value = AINSWORTH;
LASBASETOP, &choice_base_value, &flg);
if (flg != PETSC_TRUE)
if (choice_base_value == AINSWORTH)
else if (choice_base_value == DEMKOWICZ)
CHKERR simple_interface->addDomainField(
"FIELD1",
HCURL, base, 1);
CHKERR simple_interface->setFieldOrder(
"FIELD1",
order);
CHKERR simple_interface->setUp();
auto dm = simple_interface->getDM();
auto assemble_matrices_and_vectors = [&]() {
pipeline_mng->getOpDomainRhsPipeline().push_back(
pipeline_mng->getOpDomainRhsPipeline().push_back(
pipeline_mng->getOpDomainRhsPipeline().push_back(
new OpMakeHdivFromHcurl());
pipeline_mng->getOpDomainRhsPipeline().push_back(
pipeline_mng->getOpDomainRhsPipeline().push_back(
new OpAssembleVec());
pipeline_mng->getOpDomainLhsPipeline().push_back(
pipeline_mng->getOpDomainLhsPipeline().push_back(
pipeline_mng->getOpDomainLhsPipeline().push_back(
new OpMakeHdivFromHcurl());
pipeline_mng->getOpDomainLhsPipeline().push_back(
pipeline_mng->getOpDomainLhsPipeline().push_back(
new OpAssembleMat());
auto integration_rule = [](
int,
int,
int p_data) {
return 2 * p_data; };
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<MatrixDouble> ptr_values(
new MatrixDouble());
boost::shared_ptr<VectorDouble> ptr_divergence(
new VectorDouble());
pipeline_mng->getOpDomainLhsPipeline().clear();
pipeline_mng->getOpDomainRhsPipeline().clear();
pipeline_mng->getOpDomainRhsPipeline().push_back(
pipeline_mng->getOpDomainRhsPipeline().push_back(
pipeline_mng->getOpDomainRhsPipeline().push_back(
new OpMakeHdivFromHcurl());
pipeline_mng->getOpDomainRhsPipeline().push_back(
pipeline_mng->getOpDomainRhsPipeline().push_back(
pipeline_mng->getOpDomainRhsPipeline().push_back(
pipeline_mng->getOpDomainRhsPipeline().push_back(
new OpCalculateHdivVectorField<3>("FIELD1", ptr_values));
pipeline_mng->getOpDomainRhsPipeline().push_back(
new OpCalculateHdivVectorDivergence<3, 2>("FIELD1", ptr_divergence));
pipeline_mng->getOpDomainRhsPipeline().push_back(
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 .
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
@ HCURL
field with continuous tangents
#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.
int main(int argc, char *argv[])
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
OpSetInvJacHcurlFaceImpl< 2 > OpSetInvJacHcurlFace
OpCalculateJacForFaceImpl< 2 > OpCalculateJacForFace
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.
OpSetContravariantPiolaTransformFaceImpl< 2 > OpSetContravariantPiolaTransformFace
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
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.
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.