v0.13.2
Loading...
Searching...
No Matches
LinearFormsIntegrators.hpp
Go to the documentation of this file.
1/** \file LinearFormsIntegrators.hpp
2 * \brief Linear forms integrators
3 * \ingroup mofem_form
4
5*/
6
7#ifndef __LINEAR_FORMS_INTEGRATORS_HPP__
8#define __LINEAR_FORMS_INTEGRATORS_HPP__
9
10namespace MoFEM {
11
13 template <typename OpBase> struct S { S() = delete; };
15};
16
18 template <typename OpBase> struct S { S() = delete; };
20};
21
22template <int BASE_DIM, int FIELD_DIM, IntegrationType I, typename OpBase>
24
25/**
26 * @brief Integrate source
27 *
28 * @tparam OpBase
29 */
30template <typename OpBase>
32 : public OpBase {
33
34 /**
35 * @brief Construct a new Op Source Impl object
36 *
37 * @param field_name
38 * @param source_fun
39 * @param ents_ptr
40 */
41 OpSourceImpl(const std::string field_name, ScalarFun source_fun,
42 boost::shared_ptr<Range> ents_ptr = nullptr)
43 : OpBase(field_name, field_name, OpBase::OPROW, ents_ptr),
44 sourceFun(source_fun) {}
45
46 /**
47 * @brief Construct a new Op Source Impl object
48 *
49 * @param field_name
50 * @param time_fun
51 * @param source_fun
52 * @param ents_ptr
53 */
54 OpSourceImpl(const std::string field_name, TimeFun time_fun,
55 ScalarFun source_fun,
56 boost::shared_ptr<Range> ents_ptr = nullptr)
57 : OpBase(field_name, field_name, OpBase::OPROW, time_fun, ents_ptr),
58 sourceFun(source_fun) {}
59
60protected:
63};
64
65template <int FIELD_DIM, typename OpBase>
67 SourceFunctionSpecialization::S<OpBase>> : public OpBase {
68
69 /**
70 * @brief Construct a new Op Source Impl object
71 *
72 * @param field_name
73 * @param time_fun
74 * @param source_fun
75 * @param ents_ptr
76 */
77 OpSourceImpl(const std::string field_name, TimeFun time_fun,
78 VectorFun<FIELD_DIM> source_fun,
79 boost::shared_ptr<Range> ents_ptr = nullptr)
80 : OpBase(field_name, field_name, OpBase::OPROW, time_fun, ents_ptr) {}
81
82 /**
83 * @brief Construct a new Op Source Impl object
84 *
85 * @param field_name
86 * @param source_fun
87 * @param ents_ptr
88 */
89 OpSourceImpl(const std::string field_name, VectorFun<FIELD_DIM> source_fun,
90 boost::shared_ptr<Range> ents_ptr = nullptr)
91 : OpBase(field_name, field_name, OpBase::OPROW, ents_ptr),
92 sourceFun(source_fun) {}
93
94protected:
97};
98
99template <int FIELD_DIM, typename OpBase>
101 SourceFunctionSpecialization::S<OpBase>> : public OpBase {
102
103 OpSourceImpl(const std::string field_name, TimeFun time_fun,
104 VectorFun<FIELD_DIM> source_fun,
105 boost::shared_ptr<Range> ents_ptr = nullptr)
106 : OpBase(field_name, field_name, OpBase::OPROW, time_fun, ents_ptr),
107 sourceFun(source_fun) {}
108
109 OpSourceImpl(const std::string field_name, VectorFun<FIELD_DIM> source_fun,
110 boost::shared_ptr<Range> ents_ptr = nullptr)
111 : OpBase(field_name, field_name, OpBase::OPROW, ents_ptr),
112 sourceFun(source_fun) {}
113
114protected:
117};
118
119template <int BASE_DIM, int S, IntegrationType I, typename OpBase>
121
122template <int S, typename OpBase>
123struct OpBaseTimesScalarImpl<1, S, GAUSS, OpBase> : public OpBase {
124
126 const std::string field_name, boost::shared_ptr<VectorDouble> vec,
127 ScalarFun beta_coeff = [](double, double, double) constexpr { return 1; },
128 boost::shared_ptr<Range> ents_ptr = nullptr)
129 : OpBase(field_name, field_name, OpBase::OPROW, ents_ptr), sourceVec(vec),
130 betaCoeff(beta_coeff) {}
131
132protected:
134 boost::shared_ptr<VectorDouble> sourceVec;
136};
137
138template <int BASE_DIM, int FIELD_DIM, int S, IntegrationType I,
139 typename OpBase>
141
142template <int FIELD_DIM, int S, typename OpBase>
144
146 const std::string field_name, boost::shared_ptr<MatrixDouble> vec,
147 ScalarFun beta_coeff = [](double, double, double) constexpr { return 1; },
148 boost::shared_ptr<Range> ents_ptr = nullptr)
149 : OpBase(field_name, field_name, OpBase::OPROW, ents_ptr), sourceVec(vec),
150 betaCoeff(beta_coeff) {}
151
152protected:
154 boost::shared_ptr<MatrixDouble> sourceVec;
157};
158
159template <int FIELD_DIM, int S, typename OpBase>
161
163 const std::string field_name, boost::shared_ptr<MatrixDouble> vec,
164 ScalarFun beta_coeff = [](double, double, double) constexpr { return 1; },
165 boost::shared_ptr<Range> ents_ptr = nullptr)
166 : OpBase(field_name, field_name, OpBase::OPROW, ents_ptr), sourceVec(vec),
167 betaCoeff(beta_coeff) {}
168
169protected:
171 boost::shared_ptr<MatrixDouble> sourceVec;
174};
175
176template <int BASE_DIM, int FIELD_DIM, int SPACE_DIM, int S, IntegrationType I,
177 typename OpBase>
179
180template <int SPACE_DIM, int S, typename OpBase>
182 : public OpBase {
183
185
187 const std::string field_name, boost::shared_ptr<MatrixDouble> mat_vals,
188 ScalarFun beta_coeff = [](double, double, double) constexpr { return 1; },
189 boost::shared_ptr<Range> ents_ptr = nullptr)
190 : OpBase(field_name, field_name, OpBase::OPROW, ents_ptr),
191 matVals(mat_vals), betaCoeff(beta_coeff) {}
192
193protected:
194 boost::shared_ptr<MatrixDouble> matVals;
197};
198
199template <int SPACE_DIM, int S, typename OpBase>
201 : public OpBase {
202
205
207 const std::string field_name, boost::shared_ptr<MatrixDouble> mat_vals,
208 ScalarFun beta_coeff = [](double, double, double) constexpr { return 1; },
209 boost::shared_ptr<Range> ents_ptr = nullptr)
210 : OpBase(field_name, field_name, OpBase::OPROW, ents_ptr),
211 matVals(mat_vals), betaCoeff(beta_coeff) {}
212
213protected:
214 boost::shared_ptr<MatrixDouble> matVals;
217};
218
219template <int BASE_DIM, int FIELD_DIM, int SPACE_DIM, int S, IntegrationType I,
220 typename OpBase>
222
223template <int SPACE_DIM, int S, typename OpBase>
225 : public OpBase {
226
228 const std::string field_name, boost::shared_ptr<MatrixDouble> mat_vals,
229 ScalarFun beta_coeff = [](double, double, double) constexpr { return 1; })
230 : OpBase(field_name, field_name, OpBase::OPROW), matVals(mat_vals),
231 betaCoeff(beta_coeff) {}
232
233protected:
234 boost::shared_ptr<MatrixDouble> matVals;
239};
240
241template <int BASE_DIM, int FIELD_DIM, int SPACE_DIM, IntegrationType I,
242 typename OpBase>
244
245template <int FIELD_DIM, int SPACE_DIM, typename OpBase>
247 : public OpBase {
249 const std::string field_name, boost::shared_ptr<MatrixDouble> mat_vals,
250 ScalarFun beta = [](double, double, double) { return 1; },
251 boost::shared_ptr<Range> ents_ptr = nullptr)
252 : OpBase(field_name, field_name, OpBase::OPROW, ents_ptr),
253 matVals(mat_vals), betaConst(beta) {}
254
255protected:
257 boost::shared_ptr<MatrixDouble> matVals;
261};
262
263template <int SPACE_DIM, typename OpBase>
266 const std::string field_name, boost::shared_ptr<VectorDouble> vec_vals,
267 ScalarFun beta = [](double, double, double) constexpr { return 1; },
268 boost::shared_ptr<Range> ents_ptr = nullptr)
269 : OpBase(field_name, field_name, OpBase::OPROW, ents_ptr),
270 vecVals(vec_vals), betaConst(beta) {}
271
272protected:
274 boost::shared_ptr<VectorDouble> vecVals;
277};
278
279template <int FIELD_DIM, typename OpBase>
281 : public OpBase {
282
284 const std::string field_name, boost::shared_ptr<VectorDouble> vec,
285 ScalarFun beta = [](double, double, double) constexpr { return 1; },
286 boost::shared_ptr<Range> ents_ptr = nullptr)
287 : OpBase(field_name, field_name, OpBase::OPROW, ents_ptr), sourceVec(vec),
288 betaCoeff(beta) {}
289
290protected:
292 boost::shared_ptr<VectorDouble> sourceVec;
294};
295
296template <int SPACE_DIM, IntegrationType I, typename OpBase>
298
299template <int SPACE_DIM, typename OpBase>
302 boost::shared_ptr<MatrixDouble> mat_vals)
303 : OpBase(field_name, field_name, OpBase::OPROW), matVals(mat_vals) {}
304
305protected:
306 boost::shared_ptr<MatrixDouble> matVals;
309};
310
311template <int SPACE_DIM, IntegrationType I, typename OpBase>
313
314template <int SPACE_DIM, typename OpBase>
317 boost::shared_ptr<MatrixDouble> mat_vals)
318 : OpBase(field_name, field_name, OpBase::OPROW), matVals(mat_vals) {}
319
320protected:
321 boost::shared_ptr<MatrixDouble> matVals;
325};
326
327/**
328 * @brief Multiply vector times normal on the face times scalar function
329 *
330 * This operator typically will be used to evaluate natural boundary conditions
331 * for mixed formulation.
332 *
333 * @tparam BASE_DIM
334 * @tparam SPACE_DIM
335 * @tparam OpBase
336 */
337template <int SPACE_DIM, IntegrationType I, typename OpBase>
339
340/**
341 * @brief This is specialisation for sources on boundary which depends on normal
342 *
343 * @tparam BASE_DIM
344 * @tparam OpBase
345 */
346template <int FIELD_DIM, IntegrationType I, typename OpBase>
349 : public OpNormalMixVecTimesScalarImpl<FIELD_DIM, I, OpBase> {
351 OpBase>::OpNormalMixVecTimesScalarImpl;
352};
353
354template <typename OpBase>
357 const std::string field_name,
358 ScalarFun source_fun = [](double, double, double) constexpr { return 1; },
359 boost::shared_ptr<Range> ents_ptr = nullptr)
360 : OpBase(field_name, field_name, OpBase::OPROW, ents_ptr),
361 sourceFun(source_fun) {}
362
363protected:
367};
368
369template <typename OpBase>
372 const std::string field_name,
373 ScalarFun source_fun = [](double, double, double) constexpr { return 1; },
374 boost::shared_ptr<Range> ents_ptr = nullptr)
375 : OpBase(field_name, field_name, OpBase::OPROW, ents_ptr),
376 sourceFun(source_fun) {}
377
378protected:
382};
383
384template <int BASE_DIM, int FIELD_DIM, int SPACE_DIM, IntegrationType I,
385 typename OpBase>
387
388template <int SPACE_DIM, typename OpBase>
391 const std::string field_name, boost::shared_ptr<MatrixDouble> u_ptr,
392 boost::shared_ptr<MatrixDouble> y_grad_ptr,
393 ConstantFun source_fun = []() constexpr { return 1; })
394 : OpBase(field_name, field_name, OpBase::OPROW), uPtr(u_ptr),
395 yGradPtr(y_grad_ptr), alphaConstant(source_fun) {}
396
397protected:
398 boost::shared_ptr<MatrixDouble> uPtr;
399 boost::shared_ptr<MatrixDouble> yGradPtr;
402};
403
404template <int FIELD_DIM, int SPACE_DIM, typename OpBase>
406 : public OpBase {
408 const std::string field_name, boost::shared_ptr<MatrixDouble> u_ptr,
409 boost::shared_ptr<MatrixDouble> y_grad_ptr,
410 ConstantFun source_fun = []() constexpr { return 1; })
411 : OpBase(field_name, field_name, OpBase::OPROW), uPtr(u_ptr),
412 yGradPtr(y_grad_ptr), alphaConstant(source_fun) {}
413
414protected:
415 boost::shared_ptr<MatrixDouble> uPtr;
416 boost::shared_ptr<MatrixDouble> yGradPtr;
419};
420
421/**
422 * @brief Linear integrator form
423 * @ingroup mofem_forms
424 *
425 * @tparam EleOp
426 * @tparam A
427 * @tparam I
428 */
429template <typename EleOp>
430template <AssemblyType A>
431template <IntegrationType I>
432struct FormsIntegrators<EleOp>::Assembly<A>::LinearForm {
433
434 /**
435 * @brief Integrate \f$(v,f(\mathbf{x}))_\Omega\f$, f is a scalar
436 * @ingroup mofem_forms
437 *
438 * @note \f$f(\mathbf{x})\$ is scalar function of coordinates
439 *
440 * @tparam BASE_DIM
441 * @tparam FIELD_DIM
442 */
443 template <int BASE_DIM, int FIELD_DIM,
444 typename S = SourceFunctionSpecialization>
445 using OpSource =
447
448 /**
449 * @brief Vector field integrator \f$(v_i,f)_\Omega\f$, f is a vector
450 * @ingroup mofem_forms
451 *
452 * @tparam BASE_DIM
453 */
454 template <int BASE_DIM, int S = 1>
456
457 /** @deprecated use instead OpBaseTimesScalar
458 */
459 template <int BASE_DIM, int S = 1>
460 using OpBaseTimesScalarField = OpBaseTimesScalar<BASE_DIM, S>;
461
462 /**
463 * @brief Vector field integrator \f$(v,f_i)_\Omega\f$, f is a vector
464 * @ingroup mofem_forms
465 *
466 * @tparam BASE_DIM
467 * @tparam FIELD_DIM
468 * @tparam 0
469 */
470 template <int BASE_DIM, int FIELD_DIM, int S>
471 using OpBaseTimesVector =
473 //! [Source operator]
474
475 //! [Grad times tensor]
476
477 /**
478 * @brief Integrate \f$(v_{,i},f_{ij})_\Omega\f$, f is a vector
479 * @ingroup mofem_forms
480 *
481 * @note \f$f_{ij}\$ is tensor at integration points
482 *
483 * @tparam BASE_DIM
484 * @tparam FIELD_DIM
485 * @tparam SPACE_DIM
486 */
487 template <int BASE_DIM, int FIELD_DIM, int SPACE_DIM, int S = 1>
490
491 /**
492 * @brief Integrate \f$(v_{,i},f_{ij})_\Omega\f$, f is symmetric tensor
493 * @ingroup mofem_forms
494 *
495 * @note \f$f_{ij}\$ is tensor at integration points
496 *
497 * @tparam BASE_DIM
498 * @tparam FIELD_DIM
499 * @tparam SPACE_DIM
500 */
501 template <int BASE_DIM, int FIELD_DIM, int SPACE_DIM, int S = 1>
502 using OpGradTimesSymTensor =
504
505 /**
506 * @brief Integrate \f$(\lambda_{ij,j},u_{i})_\Omega\f$
507 *
508 * @tparam BASE_DIM
509 * @tparam FIELD_DIM
510 * @tparam SPACE_DIM
511 */
512 template <int BASE_DIM, int FIELD_DIM, int SPACE_DIM>
513 using OpMixDivTimesU =
515
516 /**
517 * @brief Integrate \f$(\lambda_{ij},u_{i,j})_\Omega\f$
518 *
519 * @tparam SPACE_DIM
520 */
521 template <int SPACE_DIM>
522 using OpMixTensorTimesGradU = OpMixTensorTimesGradUImpl<SPACE_DIM, I, OpBase>;
523
524 /**
525 * @brief Integrate \f$(u_{i},\lambda_{ij,j})_\Omega\f$
526 *
527 * @tparam SPACE_DIM
528 */
529 template <int SPACE_DIM>
530 using OpMixVecTimesDivLambda =
532
533 /**
534 * @brief Multiply vector times normal on the face times scalar function
535 *
536 * This operator typically will be used to evaluate natural boundary
537 * conditions for mixed formulation.
538 *
539 * @tparam BASE_DIM
540 * @tparam SPACE_DIM
541 * @tparam OpBase
542 */
543 template <int SPACE_DIM>
544 using OpNormalMixVecTimesScalar =
546
547 /**
548 * @brief Convective term
549 *
550 * \f[
551 * (v, u_i \mathbf{y}_{,i})
552 * \f]
553 * where \f$\mathbf{y}\$ can be scalar or vector.
554 *
555 * @tparam BASE_DIM
556 * @tparam FIELD_DIM
557 * @tparam SPACE_DIM
558 */
559 template <int BASE_DIM, int FIELD_DIM, int SPACE_DIM>
562};
563
564template <typename OpBase>
567 EntitiesFieldData::EntData &row_data) {
569
570 // get element volume
571 const double vol = OpBase::getMeasure();
572 // get integration weights
573 auto t_w = OpBase::getFTensor0IntegrationWeight();
574 // get base function gradient on rows
575 auto t_row_base = row_data.getFTensor0N();
576 // get coordinate at integration points
577 auto t_coords = OpBase::getFTensor1CoordsAtGaussPts();
578 // loop over integration points
579 for (int gg = 0; gg != OpBase::nbIntegrationPts; gg++) {
580 // take into account Jacobian
581 const double alpha =
582 t_w * vol * sourceFun(t_coords(0), t_coords(1), t_coords(2));
583 // loop over rows base functions
584 int rr = 0;
585 for (; rr != OpBase::nbRows; ++rr) {
586 OpBase::locF[rr] += alpha * t_row_base;
587 ++t_row_base;
588 }
589 for (; rr < OpBase::nbRowBaseFunctions; ++rr)
590 ++t_row_base;
591 ++t_coords;
592 ++t_w; // move to another integration weight
593 }
595}
596
597template <int FIELD_DIM, typename OpBase>
600 iNtegrate(EntitiesFieldData::EntData &row_data) {
603
604 // get element volume
605 const double vol = OpBase::getMeasure();
606 // get integration weights
607 auto t_w = OpBase::getFTensor0IntegrationWeight();
608 // get base function gradient on rows
609 auto t_row_base = row_data.getFTensor0N();
610 // get coordinate at integration points
611 auto t_coords = OpBase::getFTensor1CoordsAtGaussPts();
612 // loop over integration points
613 for (int gg = 0; gg != OpBase::nbIntegrationPts; gg++) {
614 // source file
615 auto t_source = sourceFun(t_coords(0), t_coords(1), t_coords(2));
616 // take into account Jacobian
617 const double alpha = t_w * vol;
618 // loop over rows base functions
619 auto t_nf = getFTensor1FromArray<FIELD_DIM, FIELD_DIM>(OpBase::locF);
620 int rr = 0;
621 for (; rr != OpBase::nbRows / FIELD_DIM; ++rr) {
622 t_nf(i) += alpha * t_row_base * t_source(i);
623 ++t_row_base;
624 ++t_nf;
625 }
626 for (; rr < OpBase::nbRowBaseFunctions; ++rr)
627 ++t_row_base;
628 ++t_coords;
629 ++t_w; // move to another integration weight
630 }
632}
633
634template <int FIELD_DIM, typename OpBase>
637 iNtegrate(EntitiesFieldData::EntData &row_data) {
640
641 const size_t nb_base_functions = row_data.getN().size2() / 3;
642 // get element volume
643 const double vol = OpBase::getMeasure();
644 // get integration weights
645 auto t_w = OpBase::getFTensor0IntegrationWeight();
646 // get base function gradient on rows
647 auto t_row_base = row_data.getFTensor1N<3>();
648 // get coordinate at integration points
649 auto t_coords = OpBase::getFTensor1CoordsAtGaussPts();
650 // loop over integration points
651 for (int gg = 0; gg != OpBase::nbIntegrationPts; gg++) {
652 // source file
653 auto t_source = sourceFun(t_coords(0), t_coords(1), t_coords(2));
654 // take into account Jacobian
655 const double alpha = t_w * vol;
656 // loop over rows base functions
657 int rr = 0;
658 for (; rr != OpBase::nbRows; ++rr) {
659 OpBase::locF[rr] += alpha * t_row_base(i) * t_source(i);
660 ++t_row_base;
661 }
662 for (; rr < nb_base_functions; ++rr)
663 ++t_row_base;
664 ++t_coords;
665 ++t_w; // move to another integration weight
666 }
668}
669
670template <int S, typename OpBase>
672 EntitiesFieldData::EntData &row_data) {
674
675 // get element volume
676 const double vol = OpBase::getMeasure();
677 // get integration weights
678 auto t_w = OpBase::getFTensor0IntegrationWeight();
679 // get base function gradient on rows
680 auto t_row_base = row_data.getFTensor0N();
681 // get vector values
682 auto t_vec = getFTensor0FromVec<S>(*sourceVec);
683
684#ifndef NDEBUG
685 if (sourceVec->size() != OpBase::nbIntegrationPts) {
686 SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
687 "Wrong number of integration points %d != %d",
688 OpBase::nbIntegrationPts, sourceVec->size());
689 }
690#endif
691
692 // get coords
693 auto t_coords = OpBase::getFTensor1CoordsAtGaussPts();
694 // loop over integration points
695 for (int gg = 0; gg != OpBase::nbIntegrationPts; gg++) {
696 // take into account Jacobian
697 const double alpha =
698 t_w * vol * betaCoeff(t_coords(0), t_coords(1), t_coords(2));
699 // loop over rows base functions
700 int rr = 0;
701 for (; rr != OpBase::nbRows; ++rr) {
702 OpBase::locF[rr] += alpha * t_row_base * t_vec;
703 ++t_row_base;
704 }
705 for (; rr < OpBase::nbRowBaseFunctions; ++rr)
706 ++t_row_base;
707 ++t_w; // move to another integration weight
708 ++t_vec;
709 ++t_coords;
710 }
712}
713
714template <int FIELD_DIM, int S, typename OpBase>
716 EntitiesFieldData::EntData &row_data) {
718
719 // get element volume
720 const double vol = OpBase::getMeasure();
721 // get integration weights
722 auto t_w = OpBase::getFTensor0IntegrationWeight();
723 // get base function gradient on rows
724 auto t_row_base = row_data.getFTensor0N();
725 // get coords
726 auto t_coords = OpBase::getFTensor1CoordsAtGaussPts();
727 // get vector values
728 auto t_vec = getFTensor1FromMat<FIELD_DIM, S>(*sourceVec);
729 // loop over integration points
730 for (int gg = 0; gg != OpBase::nbIntegrationPts; gg++) {
731 // take into account Jacobian
732 const double alpha =
733 t_w * vol * betaCoeff(t_coords(0), t_coords(1), t_coords(2));
734 // get loc vector tensor
735 auto t_nf = OpBase::template getNf<FIELD_DIM>();
736 // loop over rows base functions
737 int rr = 0;
738 for (; rr != OpBase::nbRows / FIELD_DIM; ++rr) {
739 t_nf(i) += alpha * t_row_base * t_vec(i);
740 ++t_row_base;
741 ++t_nf;
742 }
743 for (; rr < OpBase::nbRowBaseFunctions; ++rr)
744 ++t_row_base;
745 ++t_w; // move to another integration weight
746 ++t_vec;
747 ++t_coords;
748 }
750}
751
752template <int FIELD_DIM, int S, typename OpBase>
754 EntitiesFieldData::EntData &row_data) {
756
757 const size_t nb_base_functions = row_data.getN().size2() / 3;
758 // get element volume
759 const double vol = OpBase::getMeasure();
760 // get integration weights
761 auto t_w = OpBase::getFTensor0IntegrationWeight();
762 // get base function gradient on rows
763 auto t_row_base = row_data.getFTensor1N<3>();
764 // get coords
765 auto t_coords = OpBase::getFTensor1CoordsAtGaussPts();
766 // get vector values
767 auto t_vec = getFTensor1FromMat<FIELD_DIM, S>(*sourceVec);
768 // loop over integration points
769 for (int gg = 0; gg != OpBase::nbIntegrationPts; gg++) {
770 // take into account Jacobian
771 const double alpha =
772 t_w * vol * betaCoeff(t_coords(0), t_coords(1), t_coords(2));
773 // loop over rows base functions
774 int rr = 0;
775 for (; rr != OpBase::nbRows; ++rr) {
776 OpBase::locF[rr] += alpha * t_row_base(i) * t_vec(i);
777 ++t_row_base;
778 }
779 for (; rr < nb_base_functions; ++rr)
780 ++t_row_base;
781 ++t_w; // move to another integration weight
782 ++t_vec;
783 ++t_coords;
784 }
786}
787
788template <int SPACE_DIM, int S, typename OpBase>
791 EntitiesFieldData::EntData &row_data) {
793
794 // get element volume
795 const double vol = OpBase::getMeasure();
796 // get integration weights
797 auto t_w = OpBase::getFTensor0IntegrationWeight();
798 // get base function gradient on rows
799 auto t_row_grad = row_data.getFTensor1DiffN<SPACE_DIM>();
800 // get filed gradient values
801 auto t_val_grad = getFTensor1FromMat<SPACE_DIM, S>(*(matVals));
802 // get coordinate at integration points
803 auto t_coords = OpBase::getFTensor1CoordsAtGaussPts();
804 // loop over integration points
805 for (int gg = 0; gg != OpBase::nbIntegrationPts; gg++) {
806 const double beta = vol * betaCoeff(t_coords(0), t_coords(1), t_coords(2));
807 // take into account Jacobian
808 const double alpha = t_w * beta;
809 // loop over rows base functions
810 int rr = 0;
811 for (; rr != OpBase::nbRows; rr++) {
812 // calculate element of local matrix
813 OpBase::locF[rr] += alpha * (t_row_grad(i) * t_val_grad(i));
814 ++t_row_grad; // move to another element of gradient of base
815 // function on row
816 }
817 for (; rr < OpBase::nbRowBaseFunctions; ++rr)
818 ++t_row_grad;
819
820 ++t_coords;
821 ++t_val_grad;
822 ++t_w; // move to another integration weight
823 }
825}
826
827template <int SPACE_DIM, int S, typename OpBase>
830 EntitiesFieldData::EntData &row_data) {
832
833 // get element volume
834 const double vol = OpBase::getMeasure();
835 // get integration weights
836 auto t_w = OpBase::getFTensor0IntegrationWeight();
837 // get base function gradient on rows
838 auto t_row_grad = row_data.getFTensor1DiffN<SPACE_DIM>();
839 // get filed gradient values
840 auto t_val_grad = getFTensor2FromMat<SPACE_DIM, SPACE_DIM>(*(matVals));
841 // get coordinate at integration points
842 auto t_coords = OpBase::getFTensor1CoordsAtGaussPts();
843 // loop over integration points
844 for (int gg = 0; gg != OpBase::nbIntegrationPts; gg++) {
845 // take into account Jacobian
846 const double alpha =
847 t_w * vol * betaCoeff(t_coords(0), t_coords(1), t_coords(2));
848 // get rhs vector
849 auto t_nf = OpBase::template getNf<SPACE_DIM>();
850 // loop over rows base functions
851 int rr = 0;
852 for (; rr != OpBase::nbRows / SPACE_DIM; rr++) {
853 // calculate element of local matrix
854 t_nf(i) += alpha * (t_row_grad(j) * t_val_grad(i, j));
855 ++t_row_grad; // move to another element of gradient of base
856 // function on row
857 ++t_nf;
858 }
859 for (; rr < OpBase::nbRowBaseFunctions; ++rr)
860 ++t_row_grad;
861
862 ++t_coords;
863 ++t_val_grad;
864 ++t_w; // move to another integration weight
865 }
867}
868
869template <int SPACE_DIM, int S, typename OpBase>
872 EntitiesFieldData::EntData &row_data) {
874 // get element volume
875 const double vol = OpBase::getMeasure();
876 // get coords
877 auto t_coords = OpBase::getFTensor1CoordsAtGaussPts();
878 // get integration weights
879 auto t_w = OpBase::getFTensor0IntegrationWeight();
880 // get base function gradient on rows
881 auto t_row_grad = row_data.getFTensor1DiffN<SPACE_DIM>();
882 // get filed gradient values
883 auto t_val_mat = getFTensor2SymmetricFromMat<SPACE_DIM, S>(*(matVals));
884 // loop over integration points
885 for (int gg = 0; gg != OpBase::nbIntegrationPts; gg++) {
886 // take into account Jacobian
887 const double alpha =
888 t_w * vol * betaCoeff(t_coords(0), t_coords(1), t_coords(2));
889 // get rhs vector
890 auto t_nf = OpBase::template getNf<SPACE_DIM>();
891 // loop over rows base functions
892 int rr = 0;
893 for (; rr != OpBase::nbRows / SPACE_DIM; rr++) {
894 // calculate element of local matrix
895 t_nf(j) += alpha * (t_row_grad(i) * t_val_mat(i, j));
896 ++t_row_grad; // move to another element of gradient of base
897 // function on row
898 ++t_nf;
899 }
900 for (; rr < OpBase::nbRowBaseFunctions; ++rr)
901 ++t_row_grad;
902 ++t_val_mat;
903 ++t_coords;
904 ++t_w; // move to another integration weight
905 }
907}
908
909template <int FIELD_DIM, int SPACE_DIM, typename OpBase>
912 EntitiesFieldData::EntData &row_data) {
914
915 const size_t nb_base_functions = row_data.getN().size2() / 3;
916 auto t_w = this->getFTensor0IntegrationWeight();
917 // get coordinate at integration points
918 auto t_coords = OpBase::getFTensor1CoordsAtGaussPts();
919 auto t_diff_base = row_data.getFTensor2DiffN<3, SPACE_DIM>();
920 auto t_u = getFTensor1FromMat<FIELD_DIM>(*(matVals));
921
922 for (size_t gg = 0; gg != OpBase::nbIntegrationPts; ++gg) {
923
924 const double alpha = this->getMeasure() * t_w *
925 betaConst(t_coords(0), t_coords(1), t_coords(2));
926 auto t_nf = OpBase::template getNf<FIELD_DIM>();
927
928 size_t bb = 0;
929 for (; bb != this->nbRows / FIELD_DIM; ++bb) {
930 const double t_div_base = t_diff_base(j, j);
931 t_nf(i) += alpha * t_div_base * t_u(i);
932 ++t_nf;
933 ++t_diff_base;
934 }
935 for (; bb < nb_base_functions; ++bb)
936 ++t_diff_base;
937
938 ++t_u;
939 ++t_w;
940 ++t_coords;
941 }
942
944}
945
946template <int SPACE_DIM, typename OpBase>
948 EntitiesFieldData::EntData &row_data) {
950
951 const size_t nb_base_functions = row_data.getN().size2() / 3;
952 auto t_w = this->getFTensor0IntegrationWeight();
953 // get coordinate at integration points
954 auto t_coords = OpBase::getFTensor1CoordsAtGaussPts();
955 auto t_diff_base = row_data.getFTensor2DiffN<3, SPACE_DIM>();
956 auto t_u = getFTensor0FromVec(*(vecVals));
957
958 for (size_t gg = 0; gg != OpBase::nbIntegrationPts; ++gg) {
959
960 const double alpha = this->getMeasure() * t_w *
961 betaConst(t_coords(0), t_coords(1), t_coords(2));
962 ;
963
964 size_t bb = 0;
965 for (; bb != this->nbRows; ++bb) {
966 const double t_div_base = t_diff_base(j, j);
967 OpBase::locF[bb] += alpha * t_div_base * t_u;
968 ++t_diff_base;
969 }
970 for (; bb < nb_base_functions; ++bb)
971 ++t_diff_base;
972
973 ++t_u;
974 ++t_w;
975 ++t_coords;
976 }
977
979}
980
981/**
982 * @brief div U times vector
983 *
984 * \f[
985 * \delta u_j = \phi^m\delta\overline{u}^m_j\\
986 * \delta u_{j,i} = \phi^m_{,i}\delta\overline{u}^m_j\\
987 * \textrm{tr}[\delta u_{j,i}] = \delta u_{j,i}\delta_{ji}\\
988 * (\textrm{tr}[\delta u_{j,i}], v) =\\
989 * (\delta u_{j,i} \delta_{ij}, v) =\\
990 * (\delta u_{j,i}, \delta_{ij} v) =\\
991 * (\phi^m_{,i}\delta\overline{u}^m_j, \delta_{ij} v) \\
992 * f_i^m=(\phi^m_{,i}, v)
993 * \f]
994 *
995 * @tparam FIELD_DIM
996 * @tparam SPACE_DIM
997 * @tparam OpBase
998 * @param row_data
999 * @return MoFEMErrorCode
1000 */
1001template <int FIELD_DIM, typename OpBase>
1004 EntitiesFieldData::EntData &row_data) {
1007
1008 // get element volume
1009 const double vol = OpBase::getMeasure();
1010 // get integration weights
1011 auto t_w = OpBase::getFTensor0IntegrationWeight();
1012 // get coordinate at integration points
1013 auto t_coords = OpBase::getFTensor1CoordsAtGaussPts();
1014 // get base function gradient on rows
1015 auto t_row_grad = row_data.getFTensor1DiffN<FIELD_DIM>();
1016 // get vector values
1017 auto t_vec = getFTensor0FromVec(*sourceVec);
1018 // loop over integration points
1019 for (int gg = 0; gg != OpBase::nbIntegrationPts; gg++) {
1020 // take into account Jacobian
1021 const double alpha =
1022 t_w * vol * betaCoeff(t_coords(0), t_coords(1), t_coords(2));
1023 auto t_nf = getFTensor1FromArray<FIELD_DIM, FIELD_DIM>(OpBase::locF);
1024 // loop over rows base functions
1025 int rr = 0;
1026 for (; rr != OpBase::nbRows / FIELD_DIM; ++rr) {
1027 t_nf(i) += alpha * t_row_grad(i) * t_vec;
1028 ++t_row_grad;
1029 ++t_nf;
1030 }
1031 for (; rr < OpBase::nbRowBaseFunctions; ++rr)
1032 ++t_row_grad;
1033 ++t_w; // move to another integration weight
1034 ++t_vec;
1035 ++t_coords;
1036 }
1038}
1039
1040template <int SPACE_DIM, typename OpBase>
1042 EntitiesFieldData::EntData &row_data) {
1044
1045 const size_t nb_base_functions = row_data.getN().size2() / 3;
1046 auto t_w = this->getFTensor0IntegrationWeight();
1047 auto t_base = row_data.getFTensor1N<3>();
1048 auto t_grad = getFTensor2FromMat<SPACE_DIM, SPACE_DIM>(*(matVals));
1049
1050 for (size_t gg = 0; gg != OpBase::nbIntegrationPts; ++gg) {
1051
1052 const double alpha = this->getMeasure() * t_w;
1053 auto t_nf = OpBase::template getNf<SPACE_DIM>();
1054
1055 size_t bb = 0;
1056 for (; bb != this->nbRows / SPACE_DIM; ++bb) {
1057 t_nf(i) += alpha * t_base(j) * t_grad(i, j);
1058 ++t_nf;
1059 ++t_base;
1060 }
1061 for (; bb < nb_base_functions; ++bb)
1062 ++t_base;
1063
1064 ++t_grad;
1065 ++t_w;
1066 }
1067
1069}
1070
1071template <int SPACE_DIM, typename OpBase>
1073 EntitiesFieldData::EntData &row_data) {
1075
1076 const size_t nb_base_functions = row_data.getN().size2();
1077 auto t_w = this->getFTensor0IntegrationWeight();
1078 auto t_base = row_data.getFTensor0N();
1079 auto t_div = getFTensor1FromMat<SPACE_DIM>(*(matVals));
1080
1081 for (size_t gg = 0; gg != OpBase::nbIntegrationPts; ++gg) {
1082 const double alpha = this->getMeasure() * t_w;
1083 auto t_nf = OpBase::template getNf<SPACE_DIM>();
1084
1085 size_t bb = 0;
1086 for (; bb != this->nbRows / SPACE_DIM; ++bb) {
1087 t_nf(i) += alpha * t_base * t_div(i);
1088 ++t_nf;
1089 ++t_base;
1090 }
1091 for (; bb < nb_base_functions; ++bb)
1092 ++t_base;
1093
1094 ++t_div;
1095 ++t_w;
1096 }
1097
1099}
1100
1101template <typename OpBase>
1103 EntitiesFieldData::EntData &row_data) {
1105
1106 const size_t nb_base_functions = row_data.getN().size2() / 3;
1107 // get element volume
1108 const double vol = OpBase::getMeasure();
1109 // get integration weights
1110 auto t_w = OpBase::getFTensor0IntegrationWeight();
1111 // get base function gradient on rows
1112 auto t_row_base = row_data.getFTensor1N<3>();
1113 // get coordinate at integration points
1114 auto t_coords = OpBase::getFTensor1CoordsAtGaussPts();
1115 // get normal
1116 auto t_normal = OpBase::getFTensor1NormalsAtGaussPts();
1117 // loop over integration points
1118 for (int gg = 0; gg != OpBase::nbIntegrationPts; gg++) {
1119 // take into account Jacobian
1120 const double alpha = t_w * sourceFun(t_coords(0), t_coords(1), t_coords(2));
1121 // loop over rows base functions
1122 int rr = 0;
1123 for (; rr != OpBase::nbRows; ++rr) {
1124 OpBase::locF[rr] += alpha * t_row_base(i) * t_normal(i);
1125 ++t_row_base;
1126 }
1127 for (; rr < nb_base_functions; ++rr)
1128 ++t_row_base;
1129 ++t_coords;
1130 ++t_w; // move to another integration weight
1131 ++t_normal;
1132 }
1134}
1135
1136template <typename OpBase>
1138 EntitiesFieldData::EntData &row_data) {
1140
1141 const size_t nb_base_functions = row_data.getN().size2() / 3;
1142 FTensor::Tensor1<double, 3> t_z{0., 0., 1.};
1143 // get element volume
1144 const double vol = OpBase::getMeasure();
1145 // get integration weights
1146 auto t_w = OpBase::getFTensor0IntegrationWeight();
1147 // get base function gradient on rows
1148 auto t_row_base = row_data.getFTensor1N<3>();
1149 // get coordinate at integration points
1150 auto t_coords = OpBase::getFTensor1CoordsAtGaussPts();
1151 // get normal
1152 auto t_tangent = OpBase::getFTensor1TangentAtGaussPts();
1153 // loop over integration points
1154 for (int gg = 0; gg != OpBase::nbIntegrationPts; gg++) {
1155 // take into account Jacobian
1156 const double alpha = t_w * sourceFun(t_coords(0), t_coords(1), t_coords(2));
1160 t_normal(i) = FTensor::levi_civita(i, j, k) * t_tangent(j) * t_z(k);
1161 int rr = 0;
1162 for (; rr != OpBase::nbRows; ++rr) {
1163 OpBase::locF[rr] += alpha * t_row_base(i) * t_normal(i);
1164 ++t_row_base;
1165 }
1166 for (; rr < nb_base_functions; ++rr)
1167 ++t_row_base;
1168 ++t_coords;
1169 ++t_tangent;
1170 ++t_w; // move to another integration weight
1171 }
1173}
1174
1175template <int SPACE_DIM, typename OpBase>
1178 EntitiesFieldData::EntData &row_data) {
1180
1181 const size_t nb_base_functions = row_data.getN().size2();
1182 auto t_w = this->getFTensor0IntegrationWeight();
1183 auto t_base = row_data.getFTensor0N();
1184
1185 auto t_u = getFTensor1FromMat<SPACE_DIM>(*uPtr);
1186 auto t_grad_y = getFTensor1FromMat<SPACE_DIM>(*yGradPtr);
1187
1189 const double alpha_constant = alphaConstant();
1190 for (size_t gg = 0; gg != OpBase::nbIntegrationPts; ++gg) {
1191
1192 // get element volume
1193 const double vol = OpBase::getMeasure();
1194 const double c = (t_grad_y(i) * t_u(i)) * (t_w * vol * alpha_constant);
1195
1196 // get element volume
1197 int rr = 0;
1198 for (; rr != OpBase::nbRows; ++rr) {
1199 OpBase::locF[rr] += c * t_base;
1200 ++t_base;
1201 }
1202 for (; rr < OpBase::nbRowBaseFunctions; ++rr)
1203 ++t_base;
1204
1205 ++t_w; // move to another integration weight
1206 ++t_u;
1207 ++t_grad_y;
1208 }
1209
1211}
1212
1213template <int FIELD_DIM, int SPACE_DIM, typename OpBase>
1216 EntitiesFieldData::EntData &row_data) {
1218
1219 const size_t nb_base_functions = row_data.getN().size2();
1220 auto t_w = this->getFTensor0IntegrationWeight();
1221 auto t_base = row_data.getFTensor0N();
1222
1223 auto t_u = getFTensor1FromMat<SPACE_DIM>(*uPtr);
1224 auto t_grad_y = getFTensor2FromMat<FIELD_DIM, SPACE_DIM>(*yGradPtr);
1225
1228 const double alpha_constant = alphaConstant();
1229 for (size_t gg = 0; gg != OpBase::nbIntegrationPts; ++gg) {
1230
1231 // get element volume
1232 const double vol = OpBase::getMeasure();
1233
1235 t_c(J) = (t_grad_y(J, i) * t_u(i)) * (t_w * vol * alpha_constant);
1236
1237 auto t_nf = getFTensor1FromArray<FIELD_DIM, FIELD_DIM>(OpBase::locF);
1238 int rr = 0;
1239 for (; rr != OpBase::nbRows / FIELD_DIM; ++rr) {
1240 t_nf(J) += t_base * t_c(J);
1241 ++t_base;
1242 ++t_nf;
1243 }
1244 for (; rr < OpBase::nbRowBaseFunctions; ++rr)
1245 ++t_base;
1246 ++t_w; // move to another integration weight
1247
1248 ++t_u;
1249 ++t_grad_y;
1250 }
1251
1253}
1254
1255} // namespace MoFEM
1256
1257#endif // __LINEAR_FORMS_INTEGRATORS_HPP__
static Index< 'J', 3 > J
constexpr int SPACE_DIM
constexpr int FIELD_DIM
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:346
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:31
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:416
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpBaseTimesScalar< 1, 1 > OpBaseTimesScalar
IntegrationType
Form integrator integration types.
boost::function< double(const double, const double, const double)> ScalarFun
Scalar function type.
boost::function< FTensor::Tensor1< double, DIM >(const double, const double, const double)> VectorFun
Vector function type.
constexpr int BASE_DIM
FTensor::Index< 'i', SPACE_DIM > i
const double c
speed of light (cm/ns)
FTensor::Index< 'j', 3 > j
FTensor::Index< 'k', 3 > k
constexpr std::enable_if<(Dim0<=2 &&Dim1<=2), Tensor2_Expr< Levi_Civita< T >, T, Dim0, Dim1, i, j > >::type levi_civita(const Index< i, Dim0 > &, const Index< j, Dim1 > &)
levi_civita functions to make for easy adhoc use
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:10
boost::function< double(double)> TimeFun
Lambda function used to scale with time.
static auto getFTensor0FromVec(ublas::vector< T, A > &data)
Get tensor rank 0 (scalar) form data vector.
Definition: Templates.hpp:135
boost::function< double()> ConstantFun
Constant function type.
constexpr IntegrationType I
constexpr auto field_name
OpBaseImpl< PETSC, EdgeEleOp > OpBase
Definition: radiation.cpp:29
Data on single entity (This is passed as argument to DataOperator::doWork)
FTensor::Tensor2< FTensor::PackPtr< double *, Tensor_Dim0 *Tensor_Dim1 >, Tensor_Dim0, Tensor_Dim1 > getFTensor2DiffN(FieldApproximationBase base)
Get derivatives of base functions for Hdiv space.
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.
MatrixDouble & getN(const FieldApproximationBase base)
get base functions this return matrix (nb. of rows is equal to nb. of Gauss pts, nb....
FTensor::Tensor1< FTensor::PackPtr< double *, Tensor_Dim >, Tensor_Dim > getFTensor1N(FieldApproximationBase base)
Get base functions for Hdiv/Hcurl spaces.
VectorDouble locF
local entity vector
int nbRows
number of dofs on rows
int nbIntegrationPts
number of integration points
int nbRowBaseFunctions
number or row base functions
OpBaseTimesScalarImpl(const std::string field_name, boost::shared_ptr< VectorDouble > vec, ScalarFun beta_coeff=[](double, double, double) constexpr { return 1;}, boost::shared_ptr< Range > ents_ptr=nullptr)
OpBaseTimesVectorImpl(const std::string field_name, boost::shared_ptr< MatrixDouble > vec, ScalarFun beta_coeff=[](double, double, double) constexpr { return 1;}, boost::shared_ptr< Range > ents_ptr=nullptr)
OpBaseTimesVectorImpl(const std::string field_name, boost::shared_ptr< MatrixDouble > vec, ScalarFun beta_coeff=[](double, double, double) constexpr { return 1;}, boost::shared_ptr< Range > ents_ptr=nullptr)
OpConvectiveTermRhsImpl(const std::string field_name, boost::shared_ptr< MatrixDouble > u_ptr, boost::shared_ptr< MatrixDouble > y_grad_ptr, ConstantFun source_fun=[]() constexpr { return 1;})
OpConvectiveTermRhsImpl(const std::string field_name, boost::shared_ptr< MatrixDouble > u_ptr, boost::shared_ptr< MatrixDouble > y_grad_ptr, ConstantFun source_fun=[]() constexpr { return 1;})
OpGradTimesSymTensorImpl(const std::string field_name, boost::shared_ptr< MatrixDouble > mat_vals, ScalarFun beta_coeff=[](double, double, double) constexpr { return 1;})
OpGradTimesTensorImpl(const std::string field_name, boost::shared_ptr< MatrixDouble > mat_vals, ScalarFun beta_coeff=[](double, double, double) constexpr { return 1;}, boost::shared_ptr< Range > ents_ptr=nullptr)
OpGradTimesTensorImpl(const std::string field_name, boost::shared_ptr< MatrixDouble > mat_vals, ScalarFun beta_coeff=[](double, double, double) constexpr { return 1;}, boost::shared_ptr< Range > ents_ptr=nullptr)
OpMixDivTimesUImpl(const std::string field_name, boost::shared_ptr< VectorDouble > vec, ScalarFun beta=[](double, double, double) constexpr { return 1;}, boost::shared_ptr< Range > ents_ptr=nullptr)
OpMixDivTimesUImpl(const std::string field_name, boost::shared_ptr< VectorDouble > vec_vals, ScalarFun beta=[](double, double, double) constexpr { return 1;}, boost::shared_ptr< Range > ents_ptr=nullptr)
OpMixDivTimesUImpl(const std::string field_name, boost::shared_ptr< MatrixDouble > mat_vals, ScalarFun beta=[](double, double, double) { return 1;}, boost::shared_ptr< Range > ents_ptr=nullptr)
OpMixTensorTimesGradUImpl(const std::string field_name, boost::shared_ptr< MatrixDouble > mat_vals)
OpMixVecTimesDivLambdaImpl(const std::string field_name, boost::shared_ptr< MatrixDouble > mat_vals)
OpNormalMixVecTimesScalarImpl(const std::string field_name, ScalarFun source_fun=[](double, double, double) constexpr { return 1;}, boost::shared_ptr< Range > ents_ptr=nullptr)
OpNormalMixVecTimesScalarImpl(const std::string field_name, ScalarFun source_fun=[](double, double, double) constexpr { return 1;}, boost::shared_ptr< Range > ents_ptr=nullptr)
Multiply vector times normal on the face times scalar function.
OpSourceImpl(const std::string field_name, ScalarFun source_fun, boost::shared_ptr< Range > ents_ptr=nullptr)
Construct a new Op Source Impl object.
OpSourceImpl(const std::string field_name, TimeFun time_fun, ScalarFun source_fun, boost::shared_ptr< Range > ents_ptr=nullptr)
Construct a new Op Source Impl object.
OpSourceImpl(const std::string field_name, VectorFun< FIELD_DIM > source_fun, boost::shared_ptr< Range > ents_ptr=nullptr)
Construct a new Op Source Impl object.
OpSourceImpl(const std::string field_name, TimeFun time_fun, VectorFun< FIELD_DIM > source_fun, boost::shared_ptr< Range > ents_ptr=nullptr)
Construct a new Op Source Impl object.
OpSourceImpl(const std::string field_name, TimeFun time_fun, VectorFun< FIELD_DIM > source_fun, boost::shared_ptr< Range > ents_ptr=nullptr)
OpSourceImpl(const std::string field_name, VectorFun< FIELD_DIM > source_fun, boost::shared_ptr< Range > ents_ptr=nullptr)