v0.10.0
poisson_dd_H.hpp
Go to the documentation of this file.
1 #ifndef __POISSON2DTRADITIONAL_HPP__
2 #define __POISSON2DTRADITIONAL_HPP__
3 
4 #include <stdlib.h>
6 
7 #include <r_tree_testing.hpp>
8 
11 
14 
16 
18 
21 
24 
25 // common data
26 struct CommonData {
27  boost::shared_ptr<MatrixDouble> mGradPtr;
28  boost::shared_ptr<MatrixDouble> mFluxPtr;
29  boost::shared_ptr<MatrixDouble> mGradStarPtr;
30  boost::shared_ptr<MatrixDouble> mFluxStarPtr;
31 
32  boost::shared_ptr<MatrixDouble> v_field_GradPtr;
33  boost::shared_ptr<MatrixDouble> l_field_GradPtr;
34 
36  int counter;
37  bgi::rtree<point, bgi::rstar<16, 4>> rtree;
38  // RTreeDummy rTree(3);
39 };
40 
41 // const double body_source = 10;
42 typedef boost::function<double(const double, const double, const double)>
44 
45 // domain Lhs [beginning]
46 struct OpGradGradSLhs : public OpFaceEle {
47 public:
48  OpGradGradSLhs(std::string row_field_name, std::string col_field_name,
49  ScalarFunc s_p_function)
50  : OpFaceEle(row_field_name, col_field_name, OpFaceEle::OPROWCOL),
51  sFunc(s_p_function) {
52  sYmm = false;
53  }
54 
55  MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type,
56  EntityType col_type, EntData &row_data,
57  EntData &col_data) {
59 
60  const int nb_row_dofs = row_data.getIndices().size();
61  const int nb_col_dofs = col_data.getIndices().size();
62 
63  if (nb_row_dofs && nb_col_dofs) {
64 
65  locLhs.resize(nb_row_dofs, nb_col_dofs, false);
66  locLhs.clear();
67 
68  // get element area
69  const double area = getMeasure();
70 
71  // get number of integration points
72  const int nb_integration_points = getGaussPts().size2();
73  const int nb_row_base_functions = row_data.getN().size2();
74  // get integration weights
75  auto t_w = getFTensor0IntegrationWeight();
76  // get coordinates at integration point
77  auto t_coords = getFTensor1CoordsAtGaussPts();
78 
79  // get derivatives of base functions on row
80  auto t_row_diff_base = row_data.getFTensor1DiffN<2>();
81 
82  double s_term = sFunc(0, 0, 0);
83 
84  // START THE LOOP OVER INTEGRATION POINTS TO CALCULATE LOCAL MATRIX
85  for (int gg = 0; gg != nb_integration_points; gg++) {
86 
87  const double a = t_w * area;
88 
89  // get S_p (move to appropriate place after properly defined)
90  // double s_term = sFunc(t_coords(0), t_coords(1), t_coords(2));
91 
92  int rr = 0;
93  for (; rr != nb_row_dofs; ++rr) {
94  // get derivatives of base functions on column
95  auto t_col_diff_base = col_data.getFTensor1DiffN<2>(gg, 0);
96 
97  for (int cc = 0; cc != nb_col_dofs; cc++) {
98  locLhs(rr, cc) +=
99  t_row_diff_base(i) * t_col_diff_base(i) * a * s_term;
100 
101  // move to the derivatives of the next base functions on column
102  ++t_col_diff_base;
103  }
104 
105  // move to the derivatives of the next base functions on row
106  ++t_row_diff_base;
107  }
108 
109  for (; rr != nb_row_base_functions; ++rr)
110  ++t_row_diff_base;
111 
112  // move to the weight of the next integration point
113  ++t_w;
114  ++t_coords;
115  }
116 
117  // Fill value to local stiffness matrix ignoring boundary DOFs
118  CHKERR MatSetValues(getKSPB(), row_data, col_data, &locLhs(0, 0),
119  ADD_VALUES);
120  }
121 
123  }
124 
125 private:
126  boost::shared_ptr<std::vector<bool>> boundaryMarker;
129 };
130 
131 
132 // not used
133 struct OpMassSLhs : public OpFaceEle {
134 public:
135  OpMassSLhs(std::string row_field_name, std::string col_field_name,
136  ScalarFunc s_q_function)
137  : OpFaceEle(row_field_name, col_field_name, OpFaceEle::OPROWCOL),
138  sFunc(s_q_function) {
139  sYmm = false;
140  // sYmm = true;
141  }
142 
143  MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type,
144  EntityType col_type, EntData &row_data,
145  EntData &col_data) {
147 
148  const int nb_row_dofs = row_data.getIndices().size();
149  const int nb_col_dofs = col_data.getIndices().size();
150 
151  if (nb_row_dofs && nb_col_dofs) {
152 
153  locLhs.resize(nb_row_dofs, nb_col_dofs, false);
154  locLhs.clear();
155 
156  // get element area
157  const double area = getMeasure();
158 
159  // get number of integration points
160  const int nb_integration_points = getGaussPts().size2();
161  const int nb_base_functions = row_data.getN().size2();
162  // get integration weights
163  auto t_w = getFTensor0IntegrationWeight();
164  // get coordinates at integration point
165  auto t_coords = getFTensor1CoordsAtGaussPts();
166 
167  // get derivatives of base functions on row
168  auto t_row_base = row_data.getFTensor0N();
169 
170  // START THE LOOP OVER INTEGRATION POINTS TO CALCULATE LOCAL MATRIX
171  for (int gg = 0; gg != nb_integration_points; gg++) {
172 
173  const double a = t_w * area;
174 
175  // get S_p (move to appropriate place after properly defined)
176  double s_term = sFunc(t_coords(0), t_coords(1), t_coords(2));
177 
178  int rr = 0;
179  for (; rr != nb_row_dofs / 2; ++rr) {
180  // get derivatives of base functions on column
181  auto t_col_base = col_data.getFTensor0N(gg, 0);
182 
184  &locLhs(2 * rr, 0), &locLhs(2 * rr + 1, 1)};
185 
186  for (int cc = 0; cc != nb_col_dofs / 2; cc++) {
187 
188  t_a(i) += 2 * t_row_base * t_col_base * a * s_term;
189 
190  // move to the derivatives of the next base functions on column
191  ++t_col_base;
192  ++t_a;
193  }
194 
195  // move to the derivatives of the next base functions on row
196  ++t_row_base;
197  }
198 
199  for(;rr!=nb_base_functions;++rr)
200  ++t_row_base;
201 
202  // move to the weight of the next integration point
203  ++t_w;
204  }
205 
206  // FILL VALUES OF LOCAL MATRIX ENTRIES TO THE GLOBAL MATRIX
207  CHKERR MatSetValues(getKSPB(), row_data, col_data, &locLhs(0, 0),
208  ADD_VALUES);
209  }
210 
212  }
213 
214 private:
216  ScalarFunc sFunc; // add function to private part
217 };
218 
219 struct OpMGradLhs : public OpFaceEle {
220 public:
221  OpMGradLhs(std::string row_field_name, std::string col_field_name)
222  : OpFaceEle(row_field_name, col_field_name, OpFaceEle::OPROWCOL) {
223  sYmm = false;
224  }
225 
226  MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type,
227  EntityType col_type, EntData &row_data,
228  EntData &col_data) {
230 
231  const int nb_row_dofs = row_data.getIndices().size();
232  const int nb_col_dofs = col_data.getIndices().size();
233 
234  if (nb_row_dofs && nb_col_dofs) {
235 
236  const size_t nb_base_functions = row_data.getN().size2();
237  if (2 * nb_base_functions < nb_row_dofs)
238  SETERRQ(
239  PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
240  "Number of DOFs is larger than number of base functions on entity");
241 
242  locLhs.resize(nb_row_dofs, nb_col_dofs, false);
243  locLhs.clear();
244 
245  // get element area
246  const double area = getMeasure();
247 
248  // get number of integration points
249  const int nb_integration_points = getGaussPts().size2();
250  // get integration weights
251  auto t_w = getFTensor0IntegrationWeight();
252 
253  // get derivatives of base functions on row
254  auto t_row_base = row_data.getFTensor0N();
255 
256  // START THE LOOP OVER INTEGRATION POINTS TO CALCULATE LOCAL MATRIX
257  for (int gg = 0; gg != nb_integration_points; gg++) {
258 
259  const double a = t_w * area;
260 
261  int rr = 0;
262  for (; rr != nb_row_dofs / 2; ++rr) {
263 
264  // get derivatives of base functions on column
265  auto t_col_diff_base = col_data.getFTensor1DiffN<2>(gg, 0);
266 
268  &locLhs(2 * rr, 0), &locLhs(2 * rr + 1, 0)};
269 
270  for (int cc = 0; cc != nb_col_dofs; cc++) {
271 
272  t_a(i) += t_row_base * t_col_diff_base(i) * a;
273 
274  // move to the derivatives of the next base functions on column
275  ++t_col_diff_base;
276  ++t_a;
277  }
278 
279  // move to the derivatives of the next base functions on row
280  ++t_row_base;
281  }
282  for (; rr != nb_base_functions; ++rr) {
283  ++t_row_base;
284  }
285 
286  // move to the weight of the next integration point
287  ++t_w;
288  }
289 
290  // FILL VALUES OF LOCAL MATRIX ENTRIES TO THE GLOBAL MATRIX
291 
292  // Fill value to local stiffness matrix ignoring boundary DOFs
293  CHKERR MatSetValues(getKSPB(), row_data, col_data, &locLhs(0, 0),
294  ADD_VALUES);
295 
296  }
297 
299  }
300 
301 private:
302  boost::shared_ptr<std::vector<bool>> boundaryMarker;
304 };
305 
306 struct OpGradMLhs : public OpFaceEle {
307 public:
308  OpGradMLhs(std::string row_field_name, std::string col_field_name)
309  : OpFaceEle(row_field_name, col_field_name, OpFaceEle::OPROWCOL) {
310  sYmm = false;
311  }
312 
313  MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type,
314  EntityType col_type, EntData &row_data,
315  EntData &col_data) {
317 
318  const int nb_row_dofs = row_data.getIndices().size();
319  const int nb_col_dofs = col_data.getIndices().size();
320 
321  if (nb_row_dofs && nb_col_dofs) {
322 
323  const size_t nb_base_functions_row = row_data.getN().size2();
324 
325  const size_t nb_base_functions = col_data.getN().size2();
326  if (2 * nb_base_functions < nb_col_dofs)
327  SETERRQ(
328  PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
329  "Number of DOFs is larger than number of base functions on entity");
330 
331  locLhs.resize(nb_row_dofs, nb_col_dofs, false);
332  locLhs.clear();
333 
334  // get element area
335  const double area = getMeasure();
336 
337  // get number of integration points
338  const int nb_integration_points = getGaussPts().size2();
339  // get integration weights
340  auto t_w = getFTensor0IntegrationWeight();
341 
342  // get derivatives of base functions on row
343  auto t_row_diff_base = row_data.getFTensor1DiffN<2>();
344 
345 
346  // START THE LOOP OVER INTEGRATION POINTS TO CALCULATE LOCAL MATRIX
347  for (int gg = 0; gg != nb_integration_points; gg++) {
348 
349  const double a = t_w * area;
350 
351  int rr = 0;
352  for (; rr != nb_row_dofs; ++rr) {
353 
354  // get derivatives of base functions on column
355  auto t_col_base = col_data.getFTensor0N(gg, 0);
356 
358  &locLhs(rr, 0), &locLhs(rr, 1)};
359 
360  for (int cc = 0; cc != nb_col_dofs / 2; cc++) {
361 
362  t_a(i) += t_row_diff_base(i) * t_col_base * a;
363 
364  // move to the derivatives of the next base functions on column
365  ++t_col_base;
366  ++t_a;
367  }
368 
369  // move to the derivatives of the next base functions on row
370  ++t_row_diff_base;
371  }
372  for (; rr != nb_base_functions_row; ++rr) {
373  ++t_row_diff_base;
374  }
375 
376  // move to the weight of the next integration point
377  ++t_w;
378  }
379 
380  // FILL VALUES OF LOCAL MATRIX ENTRIES TO THE GLOBAL MATRIX
381  CHKERR MatSetValues(getKSPB(), row_data, col_data, &locLhs(0, 0),
382  ADD_VALUES);
383 
384  }
386  }
387 
388 private:
390 };
391 
392 // domain Lhs [end]
393 
394 // * parts of Rhs [beginning]
395 
396 struct OpDpStarRhs : public OpFaceEle {
397 public:
398  OpDpStarRhs(std::string field_name,
399  boost::shared_ptr<CommonData> common_data_ptr,
400  ScalarFunc s_p_Function)
401  : OpFaceEle(field_name, OpFaceEle::OPROW), commonDataPtr(common_data_ptr),
402  sFunc(s_p_Function) {}
403 
404  MoFEMErrorCode doWork(int side, EntityType type, EntData &data) {
406 
407  const int nb_dofs = data.getIndices().size();
408 
409  if (nb_dofs) {
410 
411  locRhs.resize(nb_dofs, false);
412  locRhs.clear();
413 
414  // get element area
415  const double area = getMeasure();
416 
417  // get number of integration points
418  const int nb_integration_points = getGaussPts().size2();
419  const int nb_base_functions = data.getN().size2();
420  // get integration weights
421  auto t_w = getFTensor0IntegrationWeight();
422  // get coordinates of the integration point
423  auto t_coords = getFTensor1CoordsAtGaussPts();
424 
425  // get base functions
426  auto t_diff_base = data.getFTensor1DiffN<2>();
427  auto s_p_gradp = sFunc(t_coords(0), t_coords(1), t_coords(2));
428  auto t_grad_star = getFTensor1FromMat<2>(*(commonDataPtr->mGradStarPtr));
429 
430  // START THE LOOP OVER INTEGRATION POINTS TO CALCULATE LOCAL VECTOR
431  for (int gg = 0; gg != nb_integration_points; ++gg) {
432 
433  const double a = t_w * area;
434 
435  int rr = 0;
436  for (; rr != nb_dofs; rr++) {
437  locRhs[rr] += t_diff_base(i) * a * s_p_gradp * t_grad_star(i);
438  // move to the next base function
439  ++t_diff_base;
440  }
441  for (; rr != nb_base_functions; ++rr)
442  ++t_diff_base;
443 
444  ++t_grad_star;
445  // move to the weight of the next integration point
446  ++t_w;
447  // move to the coordinates of the next integration point
448  ++t_coords;
449  }
450 
451  CHKERR VecSetValues(getKSPf(), data, &*locRhs.begin(), ADD_VALUES);
452 
453  }
454 
456  }
457 
458 private:
459  boost::shared_ptr<CommonData> commonDataPtr;
462 };
463 
464 struct OpqStarRhs : public OpFaceEle {
465 public:
466  OpqStarRhs(const std::string field_name,
467  boost::shared_ptr<CommonData> common_data_ptr,
468  ScalarFunc s_q_function)
469  : OpFaceEle(field_name, OpFaceEle::OPROW), sFunc(s_q_function),
470  commonDataPtr(common_data_ptr) {}
471 
472  //! [Internal force]
473  MoFEMErrorCode doWork(int side, EntityType type, EntData &data) {
475 
476  const size_t nb_dofs = data.getIndices().size();
477 
478  if (nb_dofs) {
479 
480  const size_t nb_base_functions = data.getN().size2();
481  if (2 * nb_base_functions < nb_dofs)
482  SETERRQ(
483  PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
484  "Number of DOFs is larger than number of base functions on entity");
485 
486  locRhs.resize(nb_dofs, false);
487  locRhs.clear();
488 
489  const auto nb_integration_points = getGaussPts().size2();
490  auto t_w = getFTensor0IntegrationWeight();
491  auto t_base = data.getFTensor0N();
492  auto s_q_qStar = sFunc(0, 0, 0);
493  auto t_flux_star = getFTensor1FromMat<2>(*(commonDataPtr->mFluxStarPtr));
494 
495  for (size_t gg = 0; gg != nb_integration_points; ++gg) {
496 
497  double alpha = getMeasure() * t_w;
498 
500  &locRhs[1]};
501 
502  size_t bb = 0;
503  for (; bb != nb_dofs / 2; ++bb) {
504  t_nf(i) += alpha * t_base * s_q_qStar * t_flux_star(i);
505  ++t_nf;
506  ++t_base;
507  }
508 
509  for (; bb != nb_base_functions; ++bb)
510  ++t_base;
511 
512  ++t_flux_star;
513  ++t_w;
514  }
515 
516  // FILL VALUES OF THE GLOBAL VECTOR ENTRIES FROM THE LOCAL ONES
517  CHKERR VecSetValues(getKSPf(), data, &*locRhs.begin(), ADD_VALUES);
518 
519  }
520 
522  }
523  //! [Internal force]
524 
525 private:
526  boost::shared_ptr<CommonData> commonDataPtr;
528  boost::shared_ptr<std::vector<bool>> boundaryMarker;
530 };
531 
532 // * parts of Rhs [end]
533 
534 // boundary expressions [beginning]
535 struct OpBoundaryLhs : public OpEdgeEle {
536 public:
537  OpBoundaryLhs(std::string row_field_name, std::string col_field_name)
538  : OpEdgeEle(row_field_name, col_field_name, OpEdgeEle::OPROWCOL) {
539  sYmm = false;
540  }
541 
542  MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type,
543  EntityType col_type, EntData &row_data,
544  EntData &col_data) {
546 
547  const int nb_row_dofs = row_data.getIndices().size();
548  const int nb_col_dofs = col_data.getIndices().size();
549 
550  if (nb_row_dofs && nb_col_dofs) {
551 
552  locBoundaryLhs.resize(nb_row_dofs, nb_col_dofs, false);
553  locBoundaryLhs.clear();
554 
555  // get (boundary) element length
556  const double edge = getMeasure();
557 
558  // get number of integration points
559  const int nb_integration_points = getGaussPts().size2();
560  const int nb_base_functions = row_data.getN().size2();
561  // get integration weights
562  auto t_w = getFTensor0IntegrationWeight();
563 
564  // get base functions on row
565  auto t_row_base = row_data.getFTensor0N();
566 
567  // START THE LOOP OVER INTEGRATION POINTS TO CALCULATE LOCAL MATRIX
568  for (int gg = 0; gg != nb_integration_points; gg++) {
569  const double a = t_w * edge;
570 
571  int rr = 0;
572  for (; rr != nb_row_dofs; ++rr) {
573  // get base functions on column
574  auto t_col_base = col_data.getFTensor0N(gg, 0);
575 
576  for (int cc = 0; cc != nb_col_dofs; cc++) {
577 
578  locBoundaryLhs(rr, cc) += t_row_base * t_col_base * a;
579 
580  // move to the next base functions on column
581  ++t_col_base;
582  }
583  // move to the next base functions on row
584  ++t_row_base;
585  }
586  for(;rr!=nb_base_functions;++rr)
587  ++t_row_base;
588 
589  // move to the weight of the next integration point
590  ++t_w;
591  }
592 
593  // FILL VALUES OF LOCAL MATRIX ENTRIES TO THE GLOBAL MATRIX
594  CHKERR MatSetValues(getKSPB(), row_data, col_data, &locBoundaryLhs(0, 0),
595  ADD_VALUES);
596 
597  }
598 
600  }
601 
602 private:
604 };
605 
606 struct OpBoundaryRhs : public OpEdgeEle {
607 public:
608  OpBoundaryRhs(std::string field_name, ScalarFunc boundary_function)
609  : OpEdgeEle(field_name, OpEdgeEle::OPROW),
610  boundaryFunc(boundary_function) {}
611 
612  MoFEMErrorCode doWork(int side, EntityType type, EntData &data) {
614 
615  const int nb_dofs = data.getIndices().size();
616 
617  if (nb_dofs) {
618 
619  locBoundaryRhs.resize(nb_dofs, false);
620  locBoundaryRhs.clear();
621 
622  // get (boundary) element length
623  const double edge = getMeasure();
624 
625  // get number of integration points
626  const int nb_integration_points = getGaussPts().size2();
627  const int nb_base_functions = data.getN().size2();
628  // get integration weights
629  auto t_w = getFTensor0IntegrationWeight();
630  // get coordinates at integration point
631  auto t_coords = getFTensor1CoordsAtGaussPts();
632 
633  // get base function
634  auto t_base = data.getFTensor0N();
635 
636  // START THE LOOP OVER INTEGRATION POINTS TO CALCULATE LOCAL VECTOR
637  for (int gg = 0; gg != nb_integration_points; gg++) {
638 
639  const double a = t_w * edge;
640  double boundary_term =
641  boundaryFunc(t_coords(0), t_coords(1), t_coords(2));
642 
643  int rr = 0;
644  for (; rr != nb_dofs; rr++) {
645 
646  locBoundaryRhs[rr] += t_base * boundary_term * a;
647 
648  // move to the next base function
649  ++t_base;
650  }
651  for (; rr != nb_base_functions; ++rr)
652  ++t_base;
653 
654  // move to the weight of the next integration point
655  ++t_w;
656  // move to the coordinates of the next integration point
657  ++t_coords;
658  }
659 
660  CHKERR VecSetValues(getKSPf(), data, &*locBoundaryRhs.begin(),
661  ADD_VALUES);
662 
663  }
664 
666  }
667 
668 private:
670  boost::shared_ptr<std::vector<bool>> boundaryMarker;
672 };
673 
674 // new test
675 struct OpBoundaryQbarRhs : public OpEdgeEle {
676 public:
678  std::string field_name, ScalarFunc boundary_function)
679  : OpEdgeEle(field_name, OpEdgeEle::OPROW), boundaryFunc(boundary_function) {}
680 
681  MoFEMErrorCode doWork(int side, EntityType type, EntData &data) {
683 
684  const int nb_dofs = data.getIndices().size();
685 
686  if (nb_dofs) {
687 
688  locBoundaryRhs.resize(nb_dofs, false);
689  locBoundaryRhs.clear();
690 
691  // get (boundary) element length
692  const double edge = getMeasure();
693 
694  // get number of integration points
695  const int nb_integration_points = getGaussPts().size2();
696  const int nb_base_functions = data.getN().size2();
697  // get integration weights
698  auto t_w = getFTensor0IntegrationWeight();
699  // get coordinates at integration point
700  auto t_coords = getFTensor1CoordsAtGaussPts();
701 
702  // get base function
703  auto t_base = data.getFTensor0N();
704 
705  // START THE LOOP OVER INTEGRATION POINTS TO CALCULATE LOCAL VECTOR
706  for (int gg = 0; gg != nb_integration_points; gg++) {
707 
708  const double a = t_w * edge;
709  double boundary_term =
710  boundaryFunc(t_coords(0), t_coords(1), t_coords(2));
711 
712  int rr = 0;
713  for (; rr != nb_dofs; rr++) {
714  locBoundaryRhs[rr] += t_base * boundary_term * a;
715  // move to the next base function
716  ++t_base;
717  }
718  for(;rr!=nb_base_functions;++rr)
719  ++t_base;
720 
721  // move to the weight of the next integration point
722  ++t_w;
723  // move to the coordinates of the next integration point
724  ++t_coords;
725  }
726 
727  // FILL VALUES OF LOCAL MATRIX ENTRIES TO THE GLOBAL MATRIX
728 
729  CHKERR VecSetValues(getKSPf(), data, &*locBoundaryRhs.begin(),
730  ADD_VALUES);
731  }
732 
734  }
735 
736 private:
737  boost::shared_ptr<std::vector<bool>> boundaryMarker;
740 };
741 
742 // boundary expressions [end]
743 
744 struct OpData : public OpFaceEle {
745 public:
746  OpData(std::string field_name, boost::shared_ptr<CommonData> common_data_ptr)
747  : OpFaceEle(field_name, OpFaceEle::OPROW),
748  commonDataPtr(common_data_ptr) {}
749 
750  MoFEMErrorCode doWork(int side, EntityType type, EntData &data) {
752 
753  if (type == MBVERTEX) {
754 
755  const size_t nb_gauss_pts = getGaussPts().size2();
756  commonDataPtr->mGradStarPtr->resize(2, nb_gauss_pts);
757  commonDataPtr->mFluxStarPtr->resize(2, nb_gauss_pts);
758  auto t_grad = getFTensor1FromMat<2>(*(commonDataPtr->mGradPtr));
759  auto t_flux = getFTensor1FromMat<2>(*(commonDataPtr->mFluxPtr));
760  auto t_grad_star = getFTensor1FromMat<2>(*(commonDataPtr->mGradStarPtr));
761  auto t_flux_star = getFTensor1FromMat<2>(*(commonDataPtr->mFluxStarPtr));
762 
763  for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
764 
765  t_grad_star(i) = t_grad(i);
766  t_flux_star(i) = -3 * t_grad(i);
767 
768  ++t_grad;
769  ++t_flux;
770  ++t_grad_star;
771  ++t_flux_star;
772  }
773 
774  // cerr << "AAAA" << endl;
775  // cerr << *(commonDataPtr->mGradPtr) << endl;
776  // cerr << *(commonDataPtr->mFluxPtr) << endl;
777  // cerr << *(commonDataPtr->mGradStarPtr) << endl;
778  // cerr << *(commonDataPtr->mFluxStarPtr) << endl;
779  // cerr << *(commonDataPtr->l_field_GradPtr) << endl;
780  }
782  }
783 
784 private:
785  boost::shared_ptr<CommonData> commonDataPtr;
786 };
787 
788 }; // namespace Poisson2DTraditionalDDOperators
789 
790 #endif //__POISSON2DTRADITIONAL_HPP__
OpGradGradSLhs(std::string row_field_name, std::string col_field_name, ScalarFunc s_p_function)
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
MatrixDouble & getGaussPts()
matrix of integration (Gauss) points for Volume Element
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
Operator for linear form, usually to calculate values on right hand side.
boost::shared_ptr< MatrixDouble > mGradPtr
boost::shared_ptr< std::vector< bool > > boundaryMarker
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
[Internal force]
boost::shared_ptr< MatrixDouble > mFluxPtr
OpBoundaryLhs(std::string row_field_name, std::string col_field_name)
MoFEMErrorCode VecSetValues< EssentialBcStorage >(Vec V, const DataForcesAndSourcesCore::EntData &data, const double *ptr, InsertMode iora)
Set values to vector in operator.
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, EntData &row_data, EntData &col_data)
Operator for bi-linear form, usually to calculate values on left hand side.
OpMassSLhs(std::string row_field_name, std::string col_field_name, ScalarFunc s_q_function)
ublas::matrix< double, ublas::row_major, DoubleAllocator > MatrixDouble
Definition: Types.hpp:76
EdgeElementForcesAndSourcesCoreSwitch< 0 > EdgeElementForcesAndSourcesCore
Edge finite element default.
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:485
FTensor::Tensor1< double *, Tensor_Dim > getFTensor1DiffN(const FieldApproximationBase base)
Get derivatives of base functions.
boost::shared_ptr< CommonData > commonDataPtr
[Internal force]
boost::shared_ptr< std::vector< bool > > boundaryMarker
boost::function< double(const double, const double, const double)> ScalarFunc
Edge finite elementUser is implementing own operator at Gauss points level, by own object derived fro...
MoFEMErrorCode MatSetValues< EssentialBcStorage >(Mat M, const DataForcesAndSourcesCore::EntData &row_data, const DataForcesAndSourcesCore::EntData &col_data, const double *ptr, InsertMode iora)
Set valyes to matrix in operator.
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, EntData &row_data, EntData &col_data)
Operator for bi-linear form, usually to calculate values on left hand side.
boost::shared_ptr< CommonData > commonDataPtr
FTensor::Index< 'j', 2 > j
boost::shared_ptr< CommonData > commonDataPtr
const VectorInt & getIndices() const
Get global indices of dofs on entity.
OpData(std::string field_name, boost::shared_ptr< CommonData > common_data_ptr)
bool sYmm
If true assume that matrix is symmetric structure.
boost::shared_ptr< std::vector< bool > > boundaryMarker
ForcesAndSourcesCore::UserDataOperator UserDataOperator
MoFEM::FaceElementForcesAndSourcesCore FaceEle
boost::shared_ptr< std::vector< bool > > boundaryMarker
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, EntData &row_data, EntData &col_data)
Operator for bi-linear form, usually to calculate values on left hand side.
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:67
auto getFTensor0IntegrationWeight()
Get integration weights.
DataForcesAndSourcesCore::EntData EntData
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, EntData &row_data, EntData &col_data)
Operator for bi-linear form, usually to calculate values on left hand side.
boost::shared_ptr< MatrixDouble > l_field_GradPtr
OpDpStarRhs(std::string field_name, boost::shared_ptr< CommonData > common_data_ptr, ScalarFunc s_p_Function)
boost::shared_ptr< std::vector< bool > > boundaryMarker
OpBoundaryQbarRhs(std::string field_name, ScalarFunc boundary_function)
boost::shared_ptr< MatrixDouble > mGradStarPtr
FTensor::Index< 'i', 2 > i
#define CHKERR
Inline error check.
Definition: definitions.h:604
bgi::rtree< point, bgi::rstar< 16, 4 > > rtree
OpqStarRhs(const std::string field_name, boost::shared_ptr< CommonData > common_data_ptr, ScalarFunc s_q_function)
Data on single entity (This is passed as argument to DataOperator::doWork)
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
Operator for linear form, usually to calculate values on right hand side.
boost::shared_ptr< MatrixDouble > v_field_GradPtr
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:415
ublas::vector< double, DoubleAllocator > VectorDouble
Definition: Types.hpp:74
MatrixDouble & getN(const FieldApproximationBase base)
get base functions this return matrix (nb. of rows is equal to nb. of Gauss pts, nb....
boost::shared_ptr< MatrixDouble > mFluxStarPtr
FaceElementForcesAndSourcesCoreSwitch< 0 > FaceElementForcesAndSourcesCore
Face finite element default.
OpBoundaryRhs(std::string field_name, ScalarFunc boundary_function)
OpMGradLhs(std::string row_field_name, std::string col_field_name)
FTensor::Tensor0< FTensor::PackPtr< double *, 1 > > getFTensor0N(const FieldApproximationBase base)
Get base function as Tensor0.
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, EntData &row_data, EntData &col_data)
OpGradMLhs(std::string row_field_name, std::string col_field_name)