v0.13.2
Loading...
Searching...
No Matches
minimal_surface_equation.cpp
Go to the documentation of this file.
1#include <stdlib.h>
3using namespace MoFEM;
4
5static char help[] = "...\n\n";
6
12
19
24
26
27/** \brief Integrate the domain tangent matrix (LHS)
28
29\f[
30\sum\limits_j {\left[ {\int\limits_{{\Omega _e}} {\left( {{a_n}\nabla {\phi _i}
31\cdot \nabla {\phi _j} - a_n^3\nabla {\phi _i}\left( {\nabla u \cdot \nabla
32{\phi _j}} \right)\nabla u} \right)d{\Omega _e}} } \right]\delta {U_j}} =
33\int\limits_{{\Omega _e}} {{\phi _i}fd{\Omega _e}} - \int\limits_{{\Omega _e}}
34{\nabla {\phi _i}{a_n}\nabla ud{\Omega _e}} \\
35{a_n} = \frac{1}{{{{\left( {1 +
36{{\left| {\nabla u} \right|}^2}} \right)}^{\frac{1}{2}}}}}
37\f]
38
39*/
41public:
42 OpDomainTangentMatrix(std::string row_field_name, std::string col_field_name,
43 boost::shared_ptr<MatrixDouble> field_grad_mat)
44 : AssemblyDomainEleOp(row_field_name, col_field_name,
45 DomainEleOp::OPROWCOL),
46 fieldGradMat(field_grad_mat) {}
47
48 MoFEMErrorCode iNtegrate(EntData &row_data, EntData &col_data) {
50
51 auto &locLhs = AssemblyDomainEleOp::locMat;
52
53 // get element area
54 const double area = getMeasure();
55
56 // get number of integration points
57 const int nb_integration_points = getGaussPts().size2();
58 // get integration weights
59 auto t_w = getFTensor0IntegrationWeight();
60 // get gradient of the field at integration points
61 auto t_field_grad = getFTensor1FromMat<2>(*fieldGradMat);
62
63 // get derivatives of base functions on row
64 auto t_row_diff_base = row_data.getFTensor1DiffN<2>();
65
66 // START THE LOOP OVER INTEGRATION POINTS TO CALCULATE LOCAL MATRIX
67 for (int gg = 0; gg != nb_integration_points; gg++) {
68
69 const double a = t_w * area;
70 const double an = 1. / std::sqrt(1 + t_field_grad(i) * t_field_grad(i));
71
72 for (int rr = 0; rr != AssemblyDomainEleOp::nbRows; ++rr) {
73 // get derivatives of base functions on column
74 auto t_col_diff_base = col_data.getFTensor1DiffN<2>(gg, 0);
75
76 for (int cc = 0; cc != AssemblyDomainEleOp::nbRows; cc++) {
77
78 // calculate components of the local matrix
79 locLhs(rr, cc) += (t_row_diff_base(i) * t_col_diff_base(i)) * an * a -
80 t_row_diff_base(i) *
81 (t_field_grad(i) * t_col_diff_base(i)) *
82 t_field_grad(i) * an * an * an * a;
83
84 // move to the derivatives of the next base functions on column
85 ++t_col_diff_base;
86 }
87
88 // move to the derivatives of the next base functions on row
89 ++t_row_diff_base;
90 }
91
92 // move to the weight of the next integration point
93 ++t_w;
94 // move to the gradient of the field at the next integration point
95 ++t_field_grad;
96 }
97
99 }
100
101private:
102 boost::shared_ptr<MatrixDouble> fieldGradMat;
103};
104
105/** \brief Integrate the domain residual vector (RHS)
106
107\f[
108\sum\limits_j {\left[ {\int\limits_{{\Omega _e}} {\left( {{a_n}\nabla {\phi _i}
109\cdot \nabla {\phi _j} - a_n^3\nabla {\phi _i}\left( {\nabla u \cdot \nabla
110{\phi _j}} \right)\nabla u} \right)d{\Omega _e}} } \right]\delta {U_j}} =
111\int\limits_{{\Omega _e}} {{\phi _i}fd{\Omega _e}} - \int\limits_{{\Omega _e}}
112{\nabla {\phi _i}{a_n}\nabla ud{\Omega _e}} \\
113{a_n} = \frac{1}{{{{\left( {1 +
114{{\left| {\nabla u} \right|}^2}} \right)}^{\frac{1}{2}}}}}
115\f]
116
117*/
119public:
121 boost::shared_ptr<MatrixDouble> field_grad_mat)
123 fieldGradMat(field_grad_mat) {}
124
127
128 auto &nf = AssemblyDomainEleOp::locF;
129
130 // get element area
131 const double area = getMeasure();
132
133 // get number of integration points
134 const int nb_integration_points = getGaussPts().size2();
135 // get integration weights
136 auto t_w = getFTensor0IntegrationWeight();
137 // get gradient of the field at integration points
138 auto t_field_grad = getFTensor1FromMat<2>(*fieldGradMat);
139
140 // get base functions
141 auto t_base = data.getFTensor0N();
142 // get derivatives of base functions
143 auto t_diff_base = data.getFTensor1DiffN<2>();
144
145 // START THE LOOP OVER INTEGRATION POINTS TO CALCULATE LOCAL VECTOR
146 for (int gg = 0; gg != nb_integration_points; gg++) {
147
148 const double a = t_w * area;
149 const double an = 1. / std::sqrt(1 + t_field_grad(i) * t_field_grad(i));
150
151 for (int rr = 0; rr != AssemblyDomainEleOp::nbRows; rr++) {
152
153 // calculate components of the local vector
154 // remember to use -= here due to PETSc consideration of Residual Vec
155 nf[rr] += (t_diff_base(i) * t_field_grad(i)) * an * a;
156
157 // move to the next base function
158 ++t_base;
159 // move to the derivatives of the next base function
160 ++t_diff_base;
161 }
162
163 // move to the weight of the next integration point
164 ++t_w;
165 // move to the gradient of the field at the next integration point
166 ++t_field_grad;
167 }
168
170 }
171
172private:
173 boost::shared_ptr<MatrixDouble> fieldGradMat;
174};
175
177public:
179
180 // Declaration of the main function to run analysis
182
183private:
184 // Declaration of other main functions called in runProgram()
192
193 // Function to calculate the Boundary term
194 static double boundaryFunction(const double x, const double y,
195 const double z) {
196 return sin(2 * M_PI * (x + y));
197 }
198
199 // Main interfaces
201
202 // Object to mark boundary entities for the assembling of domain elements
203 boost::shared_ptr<std::vector<unsigned char>> boundaryMarker;
205};
206
208 : mField(m_field) {}
209
212
213 readMesh();
214 setupProblem();
218 solveSystem();
220
222}
223
226
227 auto *simple = mField.getInterface<Simple>();
229 CHKERR simple->getOptions();
230 CHKERR simple->loadFile();
231
233}
234
237
238 auto *simple = mField.getInterface<Simple>();
239 CHKERR simple->addDomainField("U", H1, AINSWORTH_BERNSTEIN_BEZIER_BASE, 1);
240 CHKERR simple->addBoundaryField("U", H1, AINSWORTH_BERNSTEIN_BEZIER_BASE, 1);
241
242 int order = 3;
243 CHKERR PetscOptionsGetInt(PETSC_NULL, "", "-order", &order, PETSC_NULL);
244 CHKERR simple->setFieldOrder("U", order);
245 CHKERR simple->setUp();
246
248}
249
252
253 auto integration_rule = [](int o_row, int o_col, int approx_order) {
254 return 2 * approx_order;
255 };
256
257 auto *pipeline_mng = mField.getInterface<PipelineManager>();
258 CHKERR pipeline_mng->setDomainRhsIntegrationRule(integration_rule);
259 CHKERR pipeline_mng->setDomainLhsIntegrationRule(integration_rule);
260 CHKERR pipeline_mng->setBoundaryLhsIntegrationRule(integration_rule);
261 CHKERR pipeline_mng->setBoundaryRhsIntegrationRule(integration_rule);
262
264}
265
268
269 auto *simple = mField.getInterface<Simple>();
270
271 auto get_ents_on_mesh_skin = [&]() {
272 Range body_ents;
273 CHKERR mField.get_moab().get_entities_by_dimension(0, 2, body_ents);
274 Skinner skin(&mField.get_moab());
275 Range skin_ents;
276 CHKERR skin.find_skin(0, body_ents, false, skin_ents);
277 // filter not owned entities, those are not on boundary
278 Range boundary_ents;
279 ParallelComm *pcomm =
280 ParallelComm::get_pcomm(&mField.get_moab(), MYPCOMM_INDEX);
281 pcomm->filter_pstatus(skin_ents, PSTATUS_SHARED | PSTATUS_MULTISHARED,
282 PSTATUS_NOT, -1, &boundary_ents);
283
284 Range skin_verts;
285 mField.get_moab().get_connectivity(boundary_ents, skin_verts, true);
286 boundary_ents.merge(skin_verts);
287 boundaryEnts = boundary_ents;
288 return boundary_ents;
289 };
290
291 auto mark_boundary_dofs = [&](Range &&skin_edges) {
292 auto problem_manager = mField.getInterface<ProblemsManager>();
293 auto marker_ptr = boost::make_shared<std::vector<unsigned char>>();
294 problem_manager->markDofs(simple->getProblemName(), ROW,
295 ProblemsManager::OR, skin_edges, *marker_ptr);
296 return marker_ptr;
297 };
298
299 // Get global local vector of marked DOFs. Is global, since is set for all
300 // DOFs on processor. Is local since only DOFs on processor are in the
301 // vector. To access DOFs use local indices.
302 boundaryMarker = mark_boundary_dofs(get_ents_on_mesh_skin());
303
305}
306
309 auto add_domain_base_ops = [&](auto &pipeline) {
310 auto det_ptr = boost::make_shared<VectorDouble>();
311 auto jac_ptr = boost::make_shared<MatrixDouble>();
312 auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
313 pipeline.push_back(new OpCalculateHOJac<2>(jac_ptr));
314 pipeline.push_back(new OpInvertMatrix<2>(jac_ptr, det_ptr, inv_jac_ptr));
315 pipeline.push_back(new OpSetHOInvJacToScalarBases<2>(H1, inv_jac_ptr));
316 pipeline.push_back(new OpSetHOWeightsOnFace());
317 };
318
319 auto add_domain_lhs_ops = [&](auto &pipeline) {
320 pipeline.push_back(new OpSetBc("U", true, boundaryMarker));
321 auto grad_u_at_gauss_pts = boost::make_shared<MatrixDouble>();
322 pipeline.push_back(
323 new OpCalculateScalarFieldGradient<2>("U", grad_u_at_gauss_pts));
324 pipeline.push_back(
325 new OpDomainTangentMatrix("U", "U", grad_u_at_gauss_pts));
326 pipeline.push_back(new OpUnSetBc("U"));
327 };
328
329 auto add_domain_rhs_ops = [&](auto &pipeline) {
330 pipeline.push_back(new OpSetBc("U", true, boundaryMarker));
331 auto grad_u_at_gauss_pts = boost::make_shared<MatrixDouble>();
332 pipeline.push_back(
333 new OpCalculateScalarFieldGradient<2>("U", grad_u_at_gauss_pts));
334 pipeline.push_back(new OpDomainResidualVector("U", grad_u_at_gauss_pts));
335 pipeline.push_back(new OpUnSetBc("U"));
336 };
337
338 auto add_boundary_base_ops = [&](auto &pipeline) {};
339
340 auto add_lhs_base_ops = [&](auto &pipeline) {
341 pipeline.push_back(new OpSetBc("U", false, boundaryMarker));
342 pipeline.push_back(new OpBoundaryMass(
343 "U", "U", [](const double, const double, const double) { return 1; }));
344 pipeline.push_back(new OpUnSetBc("U"));
345 };
346 auto add_rhs_base_ops = [&](auto &pipeline) {
347 pipeline.push_back(new OpSetBc("U", false, boundaryMarker));
348 auto u_at_gauss_pts = boost::make_shared<VectorDouble>();
349 pipeline.push_back(new OpCalculateScalarFieldValues("U", u_at_gauss_pts));
350 pipeline.push_back(new OpBoundaryTimeScalarField(
351 "U", u_at_gauss_pts,
352 [](const double, const double, const double) { return 1; }));
353 pipeline.push_back(new OpBoundarySource("U", boundaryFunction));
354 pipeline.push_back(new OpUnSetBc("U"));
355 };
356
357 auto pipeline_mng = mField.getInterface<PipelineManager>();
358 add_domain_base_ops(pipeline_mng->getOpDomainLhsPipeline());
359 add_domain_base_ops(pipeline_mng->getOpDomainRhsPipeline());
360 add_domain_lhs_ops(pipeline_mng->getOpDomainLhsPipeline());
361 add_domain_rhs_ops(pipeline_mng->getOpDomainRhsPipeline());
362
363 add_boundary_base_ops(pipeline_mng->getOpBoundaryLhsPipeline());
364 add_boundary_base_ops(pipeline_mng->getOpBoundaryRhsPipeline());
365 add_lhs_base_ops(pipeline_mng->getOpBoundaryLhsPipeline());
366 add_rhs_base_ops(pipeline_mng->getOpBoundaryRhsPipeline());
367
369}
370
373
374 auto *simple = mField.getInterface<Simple>();
375
376 auto set_fieldsplit_preconditioner = [&](auto snes) {
378 KSP ksp;
379 CHKERR SNESGetKSP(snes, &ksp);
380 PC pc;
381 CHKERR KSPGetPC(ksp, &pc);
382 PetscBool is_pcfs = PETSC_FALSE;
383 PetscObjectTypeCompare((PetscObject)pc, PCFIELDSPLIT, &is_pcfs);
384
385 if (is_pcfs == PETSC_TRUE) {
386 auto bc_mng = mField.getInterface<BcManager>();
387 auto name_prb = simple->getProblemName();
388 SmartPetscObj<IS> is_all_bc;
389 CHKERR mField.getInterface<ISManager>()->isCreateProblemFieldAndRank(
390 name_prb, ROW, "U", 0, 1, is_all_bc, &boundaryEnts);
391 int is_all_bc_size;
392 CHKERR ISGetSize(is_all_bc, &is_all_bc_size);
393 MOFEM_LOG("EXAMPLE", Sev::inform)
394 << "Field split block size " << is_all_bc_size;
395 CHKERR PCFieldSplitSetIS(pc, PETSC_NULL,
396 is_all_bc); // boundary block
397 }
399 };
400
401 // Create global RHS and solution vectors
402 auto dm = simple->getDM();
403 SmartPetscObj<Vec> global_rhs, global_solution;
404 CHKERR DMCreateGlobalVector_MoFEM(dm, global_rhs);
405 global_solution = smartVectorDuplicate(global_rhs);
406
407 // Create nonlinear solver (SNES)
408 auto pipeline_mng = mField.getInterface<PipelineManager>();
409 auto solver = pipeline_mng->createSNES();
410 CHKERR SNESSetFromOptions(solver);
411 CHKERR set_fieldsplit_preconditioner(solver);
412 CHKERR SNESSetUp(solver);
413
414 // Solve the system
415 CHKERR SNESSolve(solver, global_rhs, global_solution);
416
417 // Scatter result data on the mesh
418 CHKERR DMoFEMMeshToGlobalVector(dm, global_solution, INSERT_VALUES,
419 SCATTER_REVERSE);
420
422}
423
426
427 auto post_proc = boost::make_shared<PostProcEle>(mField);
428
429 auto u_ptr = boost::make_shared<VectorDouble>();
430 post_proc->getOpPtrVector().push_back(
431 new OpCalculateScalarFieldValues("U", u_ptr));
432
434
435 post_proc->getOpPtrVector().push_back(
436
437 new OpPPMap(post_proc->getPostProcMesh(),
438 post_proc->getMapGaussPts(),
439
440 {{"U", u_ptr}},
441
442 {}, {}, {}
443
444 )
445
446 );
447
448 auto *simple = mField.getInterface<Simple>();
449 auto dm = simple->getDM();
450 CHKERR DMoFEMLoopFiniteElements(dm, simple->getDomainFEName(), post_proc);
451
452 CHKERR post_proc->writeFile("out_result.h5m");
453
455}
456
457int main(int argc, char *argv[]) {
458
459 // Initialisation of MoFEM/PETSc and MOAB data structures
460 const char param_file[] = "param_file.petsc";
462
463 // Add logging channel for example
464 auto core_log = logging::core::get();
465 core_log->add_sink(
467 LogManager::setLog("EXAMPLE");
468 MOFEM_LOG_TAG("EXAMPLE", "example")
469
470 // Error handling
471 try {
472 // Register MoFEM discrete manager in PETSc
473 DMType dm_name = "DMMOFEM";
474 CHKERR DMRegister_MoFEM(dm_name);
475
476 // Create MOAB instance
477 moab::Core mb_instance; // mesh database
478 moab::Interface &moab = mb_instance; // mesh database interface
479
480 // Create MoFEM instance
481 MoFEM::Core core(moab); // finite element database
482 MoFEM::Interface &m_field = core; // finite element interface
483
484 // Run the main analysis
485 MinimalSurfaceEqn minimal_surface_problem(m_field);
486 CHKERR minimal_surface_problem.runProgram();
487 }
489
490 // Finish work: cleaning memory, getting statistics, etc.
492
493 return 0;
494}
std::string param_file
void simple(double P1[], double P2[], double P3[], double c[], const int N)
Definition: acoustic.cpp:69
int main()
Definition: adol-c_atom.cpp:46
constexpr double a
ElementsAndOps< SPACE_DIM >::DomainEle DomainEle
ElementsAndOps< SPACE_DIM >::BoundaryEle BoundaryEle
@ ROW
Definition: definitions.h:123
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:372
@ AINSWORTH_BERNSTEIN_BEZIER_BASE
Definition: definitions.h:64
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:447
@ H1
continuous field
Definition: definitions.h:85
#define MYPCOMM_INDEX
default communicator number PCOMM
Definition: definitions.h:215
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:346
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:416
#define CHKERR
Inline error check.
Definition: definitions.h:535
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:440
auto integration_rule
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition: DMMoFEM.cpp:47
PetscErrorCode DMoFEMLoopFiniteElements(DM dm, const char fe_name[], MoFEM::FEMethod *method, CacheTupleWeakPtr cache_ptr=CacheTupleSharedPtr())
Executes FEMethod for finite elements in DM.
Definition: DMMoFEM.cpp:554
PetscErrorCode DMCreateGlobalVector_MoFEM(DM dm, Vec *g)
DMShellSetCreateGlobalVector.
Definition: DMMoFEM.cpp:1135
PetscErrorCode DMoFEMMeshToGlobalVector(DM dm, Vec g, InsertMode mode, ScatterMode scatter_mode)
set ghosted vector values on all existing mesh entities
Definition: DMMoFEM.cpp:503
SmartPetscObj< SNES > createSNES(SmartPetscObj< DM > dm=nullptr)
Create SNES (nonlinear) solver.
static LoggerType & setLog(const std::string channel)
Set ans resset chanel logger.
Definition: LogManager.cpp:364
#define MOFEM_LOG(channel, severity)
Log.
Definition: LogManager.hpp:301
#define MOFEM_LOG_TAG(channel, tag)
Tag channel.
Definition: LogManager.hpp:332
FormsIntegrators< EdgeEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpSource< 1, 1 > OpBoundarySource
Definition: helmholtz.cpp:33
FTensor::Index< 'i', 2 > i
FormsIntegrators< BoundaryEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpBaseTimesScalar< 1 > OpBoundaryTimeScalarField
static char help[]
FormsIntegrators< BoundaryEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpSource< 1, 1 > OpBoundarySource
FormsIntegrators< BoundaryEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, 1 > OpBoundaryMass
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:10
PetscErrorCode PetscOptionsGetInt(PetscOptions *, const char pre[], const char name[], PetscInt *ivalue, PetscBool *set)
SmartPetscObj< Vec > smartVectorDuplicate(SmartPetscObj< Vec > &vec)
Create duplicate vector of smart vector.
FormsIntegrators< BoundaryEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpBaseTimesScalar< 1 > OpBoundaryTimeScalarField
OpPostProcMapInMoab< SPACE_DIM, SPACE_DIM > OpPPMap
FormsIntegrators< BoundaryEleOp >::Assembly< PETSC >::BiLinearForm< G >::OpMass< 1, SPACE_DIM > OpBoundaryMass
[Only used with Hencky/nonlinear material]
Definition: plastic.cpp:71
constexpr auto field_name
static constexpr int approx_order
OpBaseImpl< PETSC, EdgeEleOp > OpBase
Definition: radiation.cpp:29
MoFEMErrorCode boundaryCondition()
MoFEMErrorCode setIntegrationRules()
static double boundaryFunction(const double x, const double y, const double z)
MinimalSurfaceEqn(MoFEM::Interface &m_field)
boost::shared_ptr< std::vector< unsigned char > > boundaryMarker
Simple interface for fast problem set-up.
Definition: BcManager.hpp:23
virtual moab::Interface & get_moab()=0
Core (interface) class.
Definition: Core.hpp:82
static MoFEMErrorCode Initialize(int *argc, char ***args, const char file[], const char help[])
Initializes the MoFEM database PETSc, MOAB and MPI.
Definition: Core.cpp:72
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition: Core.cpp:112
Deprecated interface functions.
Data on single entity (This is passed as argument to DataOperator::doWork)
FTensor::Tensor1< FTensor::PackPtr< double *, Tensor_Dim >, Tensor_Dim > getFTensor1DiffN(const FieldApproximationBase base)
Get derivatives of base functions.
FTensor::Tensor0< FTensor::PackPtr< double *, 1 > > getFTensor0N(const FieldApproximationBase base)
Get base function as Tensor0.
Section manager is used to create indexes and sections.
Definition: ISManager.hpp:23
static boost::shared_ptr< SinkType > createSink(boost::shared_ptr< std::ostream > stream_ptr, std::string comm_filter)
Create a sink object.
Definition: LogManager.cpp:279
static boost::shared_ptr< std::ostream > getStrmWorld()
Get the strm world object.
Definition: LogManager.cpp:323
VectorDouble locF
local entity vector
int nbRows
number of dofs on rows
MatrixDouble locMat
local entity block matrix
Get field gradients at integration pts for scalar filed rank 0, i.e. vector field.
Get value at integration points for scalar field.
Post post-proc data at points from hash maps.
Set indices on entities on finite element.
Set inverse jacobian to base functions.
Modify integration weights on face to take in account higher-order geometry.
PipelineManager interface.
MoFEM::FaceElementForcesAndSourcesCore FaceEle
MoFEM::EdgeElementForcesAndSourcesCore EdgeEle
Problem manager is used to build and partition problems.
MoFEMErrorCode markDofs(const std::string problem_name, RowColData rc, const enum MarkOP op, const Range ents, std::vector< unsigned char > &marker) const
Create vector with marked indices.
Simple interface for fast problem set-up.
Definition: Simple.hpp:27
intrusive_ptr for managing petsc objects
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
Integrate the domain residual vector (RHS)
boost::shared_ptr< MatrixDouble > fieldGradMat
OpDomainResidualVector(std::string field_name, boost::shared_ptr< MatrixDouble > field_grad_mat)
MoFEMErrorCode iNtegrate(EntData &data)
Class dedicated to integrate operator.
Integrate the domain tangent matrix (LHS)
boost::shared_ptr< MatrixDouble > fieldGradMat
OpDomainTangentMatrix(std::string row_field_name, std::string col_field_name, boost::shared_ptr< MatrixDouble > field_grad_mat)
MoFEMErrorCode iNtegrate(EntData &row_data, EntData &col_data)
Integrate grad-grad operator.