16static char help[] =
"...\n\n";
32 IntegrationType::GAUSS;
42template <
int FIELD_DIM>
46template <
int FIELD_DIM>
50template <
int FIELD_DIM>
54template <
int FIELD_DIM>
58template <
int FIELD_DIM>
64int main(
int argc,
char *argv[]) {
73 moab::Interface &moab = moab_core;
80 DMType dm_name =
"DMMOFEM";
84 auto core_log = logging::core::get();
117 auto post_proc = [&](
auto dm,
auto f_res,
auto out_name) {
120 boost::make_shared<PostProcBrokenMeshInMoab<DomainEle>>(m_field);
124 auto scal_vec = boost::make_shared<VectorDouble>();
125 auto vec_mat = boost::make_shared<MatrixDouble>();
127 post_proc_fe->getOpPtrVector().push_back(
129 post_proc_fe->getOpPtrVector().push_back(
133 post_proc_fe->getOpPtrVector().push_back(
137 post_proc_fe->getPostProcMesh(), post_proc_fe->getMapGaussPts(),
139 {{
"SCALAR", scal_vec}},
141 {{
"VECTOR", vec_mat}},
149 post_proc_fe->writeFile(out_name);
158 auto TestOpGradGrad = [&]() {
160 MOFEM_LOG(
"OpTester", Sev::verbose) <<
"TestOpGradGrad";
162 pip->getOpDomainLhsPipeline().clear();
163 pip->getOpDomainRhsPipeline().clear();
165 pip->setDomainLhsIntegrationRule([](
int,
int,
int o) {
return 2 *
o; });
166 pip->setDomainLhsIntegrationRule([](
int,
int,
int o) {
return 2 *
o; });
169 "SCALAR",
"SCALAR", [](
double,
double,
double) {
return 1; }));
171 "VECTOR",
"VECTOR", [](
double,
double,
double) {
return 1; }));
173 auto scl_mat = boost::make_shared<MatrixDouble>();
174 auto vec_mat = boost::make_shared<MatrixDouble>();
176 pip->getOpDomainRhsPipeline().push_back(
178 pip->getOpDomainRhsPipeline().push_back(
183 "SCALAR", scl_mat, [](
double,
double,
double) {
return 1; }));
185 "VECTOR", vec_mat, [](
double,
double,
double) {
return 1; }));
187 constexpr double eps = 1e-6;
189 auto x = opt->setRandomFields(
simple->getDM(),
190 {{
"SCALAR", {-1, 1}}, {
"VECTOR", {-1, 1}}});
191 auto diff_x = opt->setRandomFields(
192 simple->getDM(), {{
"SCALAR", {-1, 1}}, {
"VECTOR", {-1, 1}}});
194 auto diff_res = opt->checkCentralFiniteDifference(
195 simple->getDM(),
simple->getDomainFEName(), pip->getDomainRhsFE(),
202 CHKERR VecNorm(diff_res, NORM_2, &fnorm);
203 MOFEM_LOG_C(
"OpTester", Sev::inform,
"TestOpGradGrad %3.4e", fnorm);
205 constexpr double err = 1e-9;
208 "Norm of directional derivative too large");
216 opt->directionalCentralFiniteDifference(
221 "out_directional_directive_gradgrad.h5m");
228 auto TestOpConvection = [&]() {
230 MOFEM_LOG(
"OpTester", Sev::verbose) <<
"TestOpConvection";
232 pip->getOpDomainLhsPipeline().clear();
233 pip->getOpDomainRhsPipeline().clear();
235 pip->setDomainLhsIntegrationRule([](
int,
int,
int o) {
return 2 *
o; });
236 pip->setDomainLhsIntegrationRule([](
int,
int,
int o) {
return 2 *
o; });
238 auto evaluate_field = [&](
auto &pipe) {
239 auto scl_mat = boost::make_shared<MatrixDouble>();
240 auto vec_mat = boost::make_shared<MatrixDouble>();
241 auto dot_vec_vel = boost::make_shared<MatrixDouble>();
247 "VECTOR", dot_vec_vel));
248 return std::make_tuple(scl_mat, vec_mat, dot_vec_vel);
251 auto [rhs_scl_mat, rhs_vec_mat, rhs_dot_vec_vel] =
252 evaluate_field(pip->getOpDomainRhsPipeline());
253 pip->getOpDomainRhsPipeline().push_back(
255 pip->getOpDomainRhsPipeline().push_back(
259 auto [lhs_scl_mat, lhs_vec_mat, lhs_dot_vec_vel] =
260 evaluate_field(pip->getOpDomainLhsPipeline());
264 auto op_convective_term_lhs_du_scalar =
266 op_convective_term_lhs_du_scalar->feScalingFun =
267 [](
const FEMethod *fe_ptr) {
return fe_ptr->ts_a; };
268 pip->getOpDomainLhsPipeline().push_back(op_convective_term_lhs_du_scalar);
269 pip->getOpDomainLhsPipeline().push_back(
272 auto op_convective_term_lhs_du_vec =
274 op_convective_term_lhs_du_vec->feScalingFun = [](
const FEMethod *fe_ptr) {
277 pip->getOpDomainLhsPipeline().push_back(op_convective_term_lhs_du_vec);
278 pip->getOpDomainLhsPipeline().push_back(
282 constexpr double eps = 1e-6;
284 auto x = opt->setRandomFields(
simple->getDM(),
285 {{
"SCALAR", {-1, 1}}, {
"VECTOR", {-1, 1}}});
286 auto x_t = opt->setRandomFields(
287 simple->getDM(), {{
"SCALAR", {-1, 1}}, {
"VECTOR", {-1, 1}}});
288 auto diff_x = opt->setRandomFields(
289 simple->getDM(), {{
"SCALAR", {-1, 1}}, {
"VECTOR", {-1, 1}}});
291 auto diff_res = opt->checkCentralFiniteDifference(
292 simple->getDM(),
simple->getDomainFEName(), pip->getDomainRhsFE(),
299 CHKERR VecNorm(diff_res, NORM_2, &fnorm);
300 MOFEM_LOG_C(
"OpTester", Sev::inform,
"TestOpConvection %3.4e", fnorm);
302 constexpr double err = 1e-9;
305 "Norm of directional derivative too large");
314 opt->directionalCentralFiniteDifference(
316 pip->getDomainRhsFE(), x, x_t,
319 "out_directional_directive_convection.h5m");
326 CHKERR TestOpConvection();
#define MOFEM_LOG_C(channel, severity, format,...)
void simple(double P1[], double P2[], double P3[], double c[], const int N)
ElementsAndOps< SPACE_DIM >::DomainEle DomainEle
#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_ATOM_TEST_INVALID
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define CHKERR
Inline error check.
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
PetscErrorCode DMoFEMLoopFiniteElements(DM dm, const char fe_name[], MoFEM::FEMethod *method, CacheTupleWeakPtr cache_ptr=CacheTupleSharedPtr())
Executes FEMethod for finite elements in DM.
static LoggerType & setLog(const std::string channel)
Set ans resset chanel logger.
#define MOFEM_LOG(channel, severity)
Log.
#define MOFEM_LOG_TAG(channel, tag)
Tag channel.
implementation of Data Operators for Forces and Sources
FormsIntegrators< DomainEleOp >::Assembly< A >::LinearForm< I >::OpGradTimesTensor< 1, FIELD_DIM, SPACE_DIM > OpGradTimesTensor
constexpr int SPACE_DIM
[Define dimension]
FormsIntegrators< DomainEleOp >::Assembly< A >::BiLinearForm< I >::OpGradGrad< 1, FIELD_DIM, SPACE_DIM > OpGradGrad
FormsIntegrators< DomainEleOp >::Assembly< A >::BiLinearForm< I >::OpConvectiveTermLhsDy< 1, FIELD_DIM, SPACE_DIM > OpConvectiveTermLhsDy
constexpr IntegrationType I
FormsIntegrators< DomainEleOp >::Assembly< A >::LinearForm< I >::OpConvectiveTermRhs< 1, FIELD_DIM, SPACE_DIM > OpConvectiveTermRhs
FormsIntegrators< DomainEleOp >::Assembly< A >::BiLinearForm< I >::OpConvectiveTermLhsDu< 1, FIELD_DIM, SPACE_DIM > OpConvectiveTermLhsDu
OpPostProcMapInMoab< SPACE_DIM, SPACE_DIM > OpPPMap
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)
structure for User Loop Methods on finite elements
static boost::shared_ptr< SinkType > createSink(boost::shared_ptr< std::ostream > stream_ptr, std::string comm_filter)
Create a sink object.
static boost::shared_ptr< std::ostream > getStrmWorld()
Get the strm world object.
Get field gradients at integration pts for scalar filed rank 0, i.e. vector field.
Get value at integration points for scalar field.
Get field gradients at integration pts for scalar filed rank 0, i.e. vector field.
Approximate field values for given petsc vector.
Get values at integration pts for tensor filed rank 1, i.e. vector field.
Post post-proc data at points from hash maps.
Calculate directional derivative of the right hand side and compare it with tangent matrix derivative...
PipelineManager interface.
MoFEM::FaceElementForcesAndSourcesCore FaceEle
Simple interface for fast problem set-up.
intrusive_ptr for managing petsc objects
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
Volume finite element base.