v0.13.1
hcurl_check_approx_in_2d.cpp
Go to the documentation of this file.
1/**
2 * \file hcurl_check_approx_in_2d
3 * \example hcurl_check_approx_in_2d.cpp
4 *
5 * Approximate polynomial in 2D and check derivatives
6 *
7 */
8
9/* This file is part of MoFEM.
10 * MoFEM is free software: you can redistribute it and/or modify it under
11 * the terms of the GNU Lesser General Public License as published by the
12 * Free Software Foundation, either version 3 of the License, or (at your
13 * option) any later version.
14 *
15 * MoFEM is distributed in the hope that it will be useful, but WITHOUT
16 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
17 * FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
18 * License for more details.
19 *
20 * You should have received a copy of the GNU Lesser General Public
21 * License along with MoFEM. If not, see <http://www.gnu.org/licenses/>. */
22
23#include <MoFEM.hpp>
24
25using namespace MoFEM;
26
27static char help[] = "...\n\n";
28
31
32constexpr int BASE_DIM = 3;
33constexpr int SPACE_DIM = 2;
34
35/**
36 * @brief OPerator to integrate mass matrix for least square approximation
37 *
38 */
41
42/**
43 * @brief Operator to integrate the right hand side matrix for the problem
44 *
45 */
47 GAUSS>::OpSource<BASE_DIM, BASE_DIM>;
48
49constexpr double a0 = 0.0;
50constexpr double a1 = 2.0;
51constexpr double a2 = -15.0 * a0;
52constexpr double a3 = -20.0 / 6 * a1;
53constexpr double a4 = 15.0 * a0;
54constexpr double a5 = a1;
55constexpr double a6 = -a0;
57
58 static FTensor::Tensor1<double, BASE_DIM> fUn(const double x, const double y,
59 double z) {
61 // x
62 6 * a6 * std::pow(x, 5) * std::pow(y, 0) +
63 5 * a5 * std::pow(x, 4) * std::pow(y, 1) +
64 4 * a4 * std::pow(x, 3) * std::pow(y, 2) +
65 3 * a3 * std::pow(x, 2) * std::pow(y, 3) +
66 2 * a2 * std::pow(x, 1) * std::pow(y, 4) +
67 1 * a1 * std::pow(x, 0) * std::pow(y, 5),
68 // y
69 1 * a5 * std::pow(x, 5) * std::pow(y, 0) +
70 2 * a4 * std::pow(x, 4) * std::pow(y, 1) +
71 3 * a3 * std::pow(x, 3) * std::pow(y, 2) +
72 4 * a2 * std::pow(x, 2) * std::pow(y, 3) +
73 5 * a1 * std::pow(x, 1) * std::pow(y, 4) +
74 6 * a0 * std::pow(x, 0) * std::pow(y, 5),
75
76 // z
77 0.);
78 }
79
81 const double y) {
83 // x,x
84 30 * a6 * pow(x, 4) * pow(y, 0) + 20 * a5 * pow(x, 3) * pow(y, 1) +
85 12 * a4 * pow(x, 2) * pow(y, 2) + 6 * a3 * pow(x, 1) * pow(y, 3) +
86 2 * a2 * pow(x, 0) * pow(y, 4),
87 // x,y
88 5 * a5 * pow(x, 4) * pow(y, 0) + 8 * a4 * pow(x, 3) * pow(y, 1) +
89 9 * a3 * pow(x, 2) * pow(y, 2) + 8 * a2 * pow(x, 1) * pow(y, 3) +
90 5 * a1 * pow(x, 0) * pow(y, 4),
91 // y,x
92 5 * a5 * pow(x, 4) * pow(y, 0) + 8 * a4 * pow(x, 3) * pow(y, 1) +
93 9 * a3 * pow(x, 2) * pow(y, 2) + 8 * a2 * pow(x, 1) * pow(y, 3) +
94 5 * a1 * pow(x, 0) * pow(y, 4),
95 // y,y
96 2 * a4 * pow(x, 4) * pow(y, 0) + 6 * a3 * pow(x, 3) * pow(y, 1) +
97 12 * a2 * pow(x, 2) * pow(y, 2) + 20 * a1 * pow(x, 1) * pow(y, 3) +
98 30 * a0 * pow(x, 0) * pow(y, 4),
99 // z
100 0., 0.);
101 }
102
104 diff2Fun(const double x, const double y) {
106 // x,xx 0/000
107
108 30 * 4 * a6 * pow(x, 3) * pow(y, 0) +
109 20 * 3 * a5 * pow(x, 2) * pow(y, 1) +
110 12 * 2 * a4 * pow(x, 1) * pow(y, 2) +
111 6 * 1 * a3 * pow(x, 0) * pow(y, 3),
112
113 // x,xy 1/001
114
115 20 * 1 * a5 * pow(x, 3) * pow(y, 0) +
116 12 * 2 * a4 * pow(x, 2) * pow(y, 2) +
117 6 * 3 * a3 * pow(x, 1) * pow(y, 2) +
118 2 * 4 * a2 * pow(x, 0) * pow(y, 3),
119
120 // x,yx 2/010
121
122 5 * 4 * a5 * pow(x, 3) * pow(y, 0) +
123 8 * 3 * a4 * pow(x, 2) * pow(y, 1) +
124 9 * 2 * a3 * pow(x, 1) * pow(y, 2) +
125 8 * 1 * a2 * pow(x, 0) * pow(y, 3),
126
127 // x,yy 3/011
128
129 8 * 1 * a4 * pow(x, 3) * pow(y, 0) +
130 9 * 2 * a3 * pow(x, 2) * pow(y, 1) +
131 8 * 3 * a2 * pow(x, 1) * pow(y, 2) +
132 5 * 4 * a1 * pow(x, 0) * pow(y, 3),
133
134 // y,xx 4/100
135
136 5 * 4 * a5 * pow(x, 3) * pow(y, 0) +
137 8 * 3 * a4 * pow(x, 2) * pow(y, 1) +
138 9 * 2 * a3 * pow(x, 1) * pow(y, 2) +
139 8 * 1 * a2 * pow(x, 0) * pow(y, 3),
140
141 // y,xy 5/101
142
143 8 * 1 * a4 * pow(x, 3) * pow(y, 0) +
144 9 * 2 * a3 * pow(x, 2) * pow(y, 1) +
145 8 * 3 * a2 * pow(x, 1) * pow(y, 2) +
146 5 * 4 * a1 * pow(x, 0) * pow(y, 3),
147
148 // y,yx 6/110
149
150 2 * 4 * a4 * pow(x, 3) * pow(y, 0) +
151 6 * 3 * a3 * pow(x, 2) * pow(y, 1) +
152 12 * 2 * a2 * pow(x, 1) * pow(y, 2) +
153 20 * 1 * a1 * pow(x, 0) * pow(y, 3),
154
155 // y,yy 7/111
156
157 6 * 1 * a3 * pow(x, 3) * pow(y, 0) +
158 12 * 2 * a2 * pow(x, 2) * pow(y, 1) +
159 20 * 3 * a1 * pow(x, 1) * pow(y, 2) +
160 30 * 4 * a0 * pow(x, 0) * pow(y, 3),
161
162 // z,xx 8/200
163 0.,
164 // z,xy 9/201
165 0.,
166 // z,yx 10/210
167 0.,
168 // z,yy 11/211
169 0.);
170 }
171};
172
174 boost::shared_ptr<MatrixDouble> ptrVals;
175 boost::shared_ptr<VectorDouble> ptrDiv;
176 boost::shared_ptr<MatrixDouble> ptrGrad;
177 boost::shared_ptr<MatrixDouble> ptrHess;
178
179 OpCheckValsDiffVals(boost::shared_ptr<MatrixDouble> ptr_vals,
180 boost::shared_ptr<VectorDouble> ptr_div,
181 boost::shared_ptr<MatrixDouble> ptr_grad,
182 boost::shared_ptr<MatrixDouble> ptr_hess)
183 : FaceEleOp("FIELD1", OPROW), ptrVals(ptr_vals), ptrDiv(ptr_div),
184 ptrGrad(ptr_grad), ptrHess(ptr_hess) {}
185
189
193 const double eps = 1e-6;
194 if (type == MBEDGE && side == 0) {
195 const int nb_gauss_pts = data.getN().size1();
196
197 auto t_vals_from_op = getFTensor1FromMat<3>(*ptrVals);
198 auto t_div_from_op = getFTensor0FromVec(*ptrDiv);
199 auto t_grad_from_op = getFTensor2FromMat<3, 2>(*ptrGrad);
200 auto t_hess_from_op = getFTensor3FromMat<3, 2, 2>(*ptrHess);
201
202 for (int gg = 0; gg != nb_gauss_pts; gg++) {
203 const double x = getCoordsAtGaussPts()(gg, 0);
204 const double y = getCoordsAtGaussPts()(gg, 1);
205
206 // Check approximation
208 delta_val(i) = t_vals_from_op(i) - ApproxFunctions::fUn(x, y, 0)(i);
209 double err_val = sqrt(delta_val(i) * delta_val(i));
210 if (err_val > eps)
211 SETERRQ1(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
212 "Wrong value %4.3e", err_val);
213
214 FTensor::Tensor2<double, 3, 2> delta_diff_val;
215 delta_diff_val(i, j) =
216 t_grad_from_op(i, j) - ApproxFunctions::diffFun(x, y)(i, j);
217 double err_diff_val = sqrt(delta_diff_val(i, j) * delta_diff_val(i, j));
218 if (err_diff_val > eps)
219 SETERRQ1(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
220 "Wrong derivative of value %4.3e", err_diff_val);
221
222 double div = t_grad_from_op(0, 0) + t_grad_from_op(1, 1);
223 double err_div = div - t_div_from_op;
224 if (err_div > eps)
225 SETERRQ3(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
226 "Wrong divergence from operator %4.3e (%4.3e != %4.3e)",
227 err_div, div, t_div_from_op);
228
229 FTensor::Tensor3<double, 3, 2, 2> delta_diff2_val;
230 delta_diff2_val(i, j, k) =
231 t_hess_from_op(i, j, k) - ApproxFunctions::diff2Fun(x, y)(i, j, k);
232 double hess_diff_error =
233 sqrt(delta_diff2_val(i, j, k) * delta_diff2_val(i, j, k));
234 if (hess_diff_error > eps)
235 SETERRQ1(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
236 "Wrong hessian from operator %4.3e", hess_diff_error);
237
238 ++t_vals_from_op;
239 ++t_div_from_op;
240 ++t_grad_from_op;
241 ++t_hess_from_op;
242 }
243 }
245 }
246};
247
248int main(int argc, char *argv[]) {
249
250 MoFEM::Core::Initialize(&argc, &argv, (char *)0, help);
251
252 try {
253
254 DMType dm_name = "DMMOFEM";
255 CHKERR DMRegister_MoFEM(dm_name);
256
257 moab::Core mb_instance;
258 moab::Interface &moab = mb_instance;
259
260 // Create MoFEM instance
261 MoFEM::Core core(moab);
262 MoFEM::Interface &m_field = core;
263
264 Simple *simple_interface = m_field.getInterface<Simple>();
265 PipelineManager *pipeline_mng = m_field.getInterface<PipelineManager>();
266 CHKERR simple_interface->getOptions();
267 CHKERR simple_interface->loadFile("", "rectangle_tri.h5m");
268
269 // Declare elements
270 enum bases { AINSWORTH, DEMKOWICZ, LASBASETOP };
271 const char *list_bases[] = {"ainsworth", "demkowicz"};
272 PetscBool flg;
273 PetscInt choice_base_value = AINSWORTH;
274 CHKERR PetscOptionsGetEList(PETSC_NULL, NULL, "-base", list_bases,
275 LASBASETOP, &choice_base_value, &flg);
276
277 if (flg != PETSC_TRUE)
278 SETERRQ(PETSC_COMM_SELF, MOFEM_IMPOSIBLE_CASE, "base not set");
280 if (choice_base_value == AINSWORTH)
282 else if (choice_base_value == DEMKOWICZ)
284
285 CHKERR simple_interface->addDomainField("FIELD1", HCURL, base, 1);
286 constexpr int order = 5;
287 CHKERR simple_interface->setFieldOrder("FIELD1", order);
288 CHKERR simple_interface->setUp();
289 auto dm = simple_interface->getDM();
290
291 MatrixDouble vals, diff_vals;
292
293 auto assemble_matrices_and_vectors = [&]() {
295 auto jac_ptr = boost::make_shared<MatrixDouble>();
296 auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
297 auto det_ptr = boost::make_shared<VectorDouble>();
298
299 pipeline_mng->getOpDomainRhsPipeline().push_back(
300 new OpCalculateHOJac<SPACE_DIM>(jac_ptr));
301 pipeline_mng->getOpDomainRhsPipeline().push_back(
302 new OpInvertMatrix<SPACE_DIM>(jac_ptr, det_ptr, inv_jac_ptr));
303 pipeline_mng->getOpDomainRhsPipeline().push_back(
304 new OpMakeHdivFromHcurl());
305 pipeline_mng->getOpDomainRhsPipeline().push_back(
307 pipeline_mng->getOpDomainRhsPipeline().push_back(
308 new OpDomainSource("FIELD1", ApproxFunctions::fUn));
309
310 pipeline_mng->getOpDomainLhsPipeline().push_back(
311 new OpCalculateHOJac<SPACE_DIM>(jac_ptr));
312 pipeline_mng->getOpDomainLhsPipeline().push_back(
313 new OpInvertMatrix<SPACE_DIM>(jac_ptr, det_ptr, inv_jac_ptr));
314 pipeline_mng->getOpDomainLhsPipeline().push_back(
315 new OpMakeHdivFromHcurl());
316 pipeline_mng->getOpDomainLhsPipeline().push_back(
318 pipeline_mng->getOpDomainLhsPipeline().push_back(
319
320 new OpDomainMass("FIELD1", "FIELD1",
321 [](double, double, double) { return 1; })
322
323 );
324
325 auto integration_rule = [](int, int, int p_data) { return 2 * p_data; };
326 CHKERR pipeline_mng->setDomainRhsIntegrationRule(integration_rule);
327 CHKERR pipeline_mng->setDomainLhsIntegrationRule(integration_rule);
328
330 };
331
332 auto solve_problem = [&] {
334 auto solver = pipeline_mng->createKSP();
335 CHKERR KSPSetFromOptions(solver);
336 CHKERR KSPSetUp(solver);
337
338 auto D = smartCreateDMVector(dm);
339 auto F = smartVectorDuplicate(D);
340
341 CHKERR KSPSolve(solver, F, D);
342 CHKERR VecGhostUpdateBegin(D, INSERT_VALUES, SCATTER_FORWARD);
343 CHKERR VecGhostUpdateEnd(D, INSERT_VALUES, SCATTER_FORWARD);
344 CHKERR DMoFEMMeshToLocalVector(dm, D, INSERT_VALUES, SCATTER_REVERSE);
346 };
347
348 auto check_solution = [&] {
350
351 pipeline_mng->getOpDomainLhsPipeline().clear();
352 pipeline_mng->getOpDomainRhsPipeline().clear();
353
354 auto ptr_values = boost::make_shared<MatrixDouble>();
355 auto ptr_divergence = boost::make_shared<VectorDouble>();
356 auto ptr_grad = boost::make_shared<MatrixDouble>();
357 auto ptr_hessian = boost::make_shared<MatrixDouble>();
358
359 auto jac_ptr = boost::make_shared<MatrixDouble>();
360 auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
361 auto det_ptr = boost::make_shared<VectorDouble>();
362
363 // Change H-curl to H-div in 2D, and apply Piola transform
364 pipeline_mng->getOpDomainRhsPipeline().push_back(
365 new OpCalculateHOJac<SPACE_DIM>(jac_ptr));
366 pipeline_mng->getOpDomainRhsPipeline().push_back(
367 new OpInvertMatrix<SPACE_DIM>(jac_ptr, det_ptr, inv_jac_ptr));
368 pipeline_mng->getOpDomainRhsPipeline().push_back(
369 new OpMakeHdivFromHcurl());
370 pipeline_mng->getOpDomainRhsPipeline().push_back(
372 pipeline_mng->getOpDomainRhsPipeline().push_back(
373 new OpSetInvJacHcurlFace(inv_jac_ptr));
374
375 // Evaluate base function second derivative
376 auto base_mass = boost::make_shared<MatrixDouble>();
377 auto data_l2 = boost::make_shared<EntitiesFieldData>(MBENTITYSET);
378
379 pipeline_mng->getOpDomainRhsPipeline().push_back(
380 new OpBaseDerivativesMass<BASE_DIM>(base_mass, data_l2, base, L2));
381 pipeline_mng->getOpDomainRhsPipeline().push_back(
383 inv_jac_ptr));
384 pipeline_mng->getOpDomainRhsPipeline().push_back(
385 new OpBaseDerivativesNext<BASE_DIM>(BaseDerivatives::SecondDerivative,
386 base_mass, data_l2, base, HCURL));
387
388 // Calculate field values at integration points
389 pipeline_mng->getOpDomainRhsPipeline().push_back(
390 new OpCalculateHVecVectorField<BASE_DIM>("FIELD1", ptr_values));
391 // Gradient
392 pipeline_mng->getOpDomainRhsPipeline().push_back(
394 ptr_grad));
395 // Hessian
396 pipeline_mng->getOpDomainRhsPipeline().push_back(
398 "FIELD1", ptr_divergence));
399 pipeline_mng->getOpDomainRhsPipeline().push_back(
401 ptr_hessian));
402
403 pipeline_mng->getOpDomainRhsPipeline().push_back(new OpCheckValsDiffVals(
404 ptr_values, ptr_divergence, ptr_grad, ptr_hessian));
405
406 CHKERR pipeline_mng->loopFiniteElements();
407
409 };
410
411 CHKERR assemble_matrices_and_vectors();
412 CHKERR solve_problem();
413 CHKERR check_solution();
414 }
416
418}
static const double eps
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpSource< 1, FIELD_DIM > OpDomainSource
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, FIELD_DIM > OpDomainMass
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:385
FieldApproximationBase
approximation base
Definition: definitions.h:71
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:73
@ DEMKOWICZ_JACOBI_BASE
Definition: definitions.h:79
@ L2
field with C-1 continuity
Definition: definitions.h:101
@ HCURL
field with continuous tangents
Definition: definitions.h:99
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:359
@ MOFEM_IMPOSIBLE_CASE
Definition: definitions.h:48
@ MOFEM_ATOM_TEST_INVALID
Definition: definitions.h:53
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:429
#define CHKERR
Inline error check.
Definition: definitions.h:548
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, SPACE_DIM > OpMass
auto integration_rule
PetscErrorCode DMoFEMMeshToLocalVector(DM dm, Vec l, InsertMode mode, ScatterMode scatter_mode)
set local (or ghosted) vector values on mesh for partition only
Definition: DMMMoFEM.cpp:482
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition: DMMMoFEM.cpp:59
auto smartCreateDMVector(DM dm)
Get smart vector from DM.
Definition: DMMoFEM.hpp:977
constexpr double a5
constexpr double a2
int main(int argc, char *argv[])
FormsIntegrators< FaceEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpSource< BASE_DIM, BASE_DIM > OpDomainSource
Operator to integrate the right hand side matrix for the problem.
static char help[]
constexpr double a4
constexpr int SPACE_DIM
constexpr double a0
constexpr int BASE_DIM
constexpr double a3
constexpr double a1
constexpr double a6
FormsIntegrators< FaceEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< BASE_DIM, BASE_DIM > OpDomainMass
OPerator to integrate mass matrix for least square approximation.
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:67
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:21
OpSetInvJacHcurlFaceImpl< 2 > OpSetInvJacHcurlFace
SmartPetscObj< Vec > smartVectorDuplicate(SmartPetscObj< Vec > &vec)
Create duplicate vector of smart vector.
CoreTmp< 0 > Core
Definition: Core.hpp:1096
OpSetContravariantPiolaTransformOnFace2DImpl< 2 > OpSetContravariantPiolaTransformOnFace2D
static auto getFTensor0FromVec(ublas::vector< T, A > &data)
Get tensor rank 0 (scalar) form data vector.
Definition: Templates.hpp:149
PetscErrorCode PetscOptionsGetEList(PetscOptions *, const char pre[], const char name[], const char *const *list, PetscInt next, PetscInt *value, PetscBool *set)
DeprecatedCoreInterface Interface
Definition: Interface.hpp:1965
const double D
diffusivity
static FTensor::Tensor2< double, BASE_DIM, SPACE_DIM > diffFun(const double x, const double y)
static FTensor::Tensor3< double, BASE_DIM, SPACE_DIM, SPACE_DIM > diff2Fun(const double x, const double y)
static FTensor::Tensor1< double, BASE_DIM > fUn(const double x, const double y, double z)
Core (interface) class.
Definition: Core.hpp:92
static MoFEMErrorCode Initialize(int *argc, char ***args, const char file[], const char help[])
Initializes the MoFEM database PETSc, MOAB and MPI.
Definition: Core.cpp:85
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition: Core.cpp:125
Deprecated interface functions.
Data on single entity (This is passed as argument to DataOperator::doWork)
MatrixDouble & getN(const FieldApproximationBase base)
get base functions this return matrix (nb. of rows is equal to nb. of Gauss pts, nb....
MatrixDouble & getCoordsAtGaussPts()
Gauss points and weight, matrix (nb. of points x 3)
@ OPROW
operator doWork function is executed on FE rows
Get vector field for H-div approximation.
Calculate gradient of vector field.
Calculate gradient of vector field.
Calculate divergence of vector field.
Make Hdiv space from Hcurl space in 2d.
PipelineManager interface.
Simple interface for fast problem set-up.
Definition: Simple.hpp:35
MoFEMErrorCode addDomainField(const std::string &name, const FieldSpace space, const FieldApproximationBase base, const FieldCoefficientsNumber nb_of_coefficients, const TagType tag_type=MB_TAG_SPARSE, const enum MoFEMTypes bh=MF_ZERO, int verb=-1)
Add field on domain.
Definition: Simple.cpp:378
MoFEMErrorCode getOptions()
get options
Definition: Simple.cpp:294
MoFEMErrorCode getDM(DM *dm)
Get DM.
Definition: Simple.cpp:810
MoFEMErrorCode loadFile(const std::string options, const std::string mesh_file_name)
Load mesh file.
Definition: Simple.cpp:308
MoFEMErrorCode setFieldOrder(const std::string field_name, const int order, const Range *ents=NULL)
Set field order.
Definition: Simple.cpp:593
MoFEMErrorCode setUp(const PetscBool is_partitioned=PETSC_TRUE)
Setup problem.
Definition: Simple.cpp:753
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
FTensor::Index< 'i', 3 > i
boost::shared_ptr< MatrixDouble > ptrHess
boost::shared_ptr< MatrixDouble > ptrGrad
boost::shared_ptr< VectorDouble > ptrDiv
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
OpCheckValsDiffVals(boost::shared_ptr< MatrixDouble > ptr_vals, boost::shared_ptr< VectorDouble > ptr_div, boost::shared_ptr< MatrixDouble > ptr_grad, boost::shared_ptr< MatrixDouble > ptr_hess)
FTensor::Index< 'j', 2 > j
boost::shared_ptr< MatrixDouble > ptrVals
FTensor::Index< 'k', 2 > k