12 #ifndef MFEM_TEMPLATE_BILININTEG
13 #define MFEM_TEMPLATE_BILININTEG
15 #include "../config/tconfig.hpp"
26 template <
typename coeff_t,
template<
int,
int,
typename>
class kernel_t>
32 template <
int SDim,
int Dim,
typename complex_t>
33 struct kernel {
typedef kernel_t<SDim,Dim,complex_t>
type; };
43 template <
int SDim,
int Dim,
typename complex_t>
68 template <
typename IR,
typename coeff_t,
int NE>
81 template <
typename T_result_t,
typename Q_t,
typename q_t,
83 static inline MFEM_ALWAYS_INLINE
84 void Action(
const int k,
const T_result_t &F,
85 const Q_t &Q,
const q_t &q, S_data_t &R)
87 typedef typename T_result_t::Jt_type::data_type
real_t;
88 const int M = S_data_t::eval_type::qpts;
89 const int NC = S_data_t::eval_type::vdim;
90 MFEM_STATIC_ASSERT(T_result_t::Jt_type::layout_type::dim_1 == M,
91 "incompatible dimensions");
92 MFEM_FLOPS_ADD(M*(1+NC));
93 for (
int i = 0; i < M; i++)
96 Q.get(q,i,k) * TDet<real_t>(F.Jt.layout.ind14(i,k), F.Jt);
97 for (
int j = 0; j < NC; j++)
99 R.val_qpts(i,j,k) *= wi;
111 template <
typename T_result_t,
typename Q_t,
typename q_t,
int qpts>
112 static inline MFEM_ALWAYS_INLINE
116 typedef typename T_result_t::Jt_type::data_type
real_t;
118 const int M = T_result_t::Jt_type::layout_type::dim_1;
119 MFEM_STATIC_ASSERT(qpts == M,
"incompatible dimensions");
121 for (
int i = 0; i < M; i++)
123 A[i] = Q.get(q,i,k) * TDet<real_t>(F.Jt.layout.ind14(i,k), F.Jt);
132 template <
int qpts,
typename S_data_t>
133 static inline MFEM_ALWAYS_INLINE
136 const int M = S_data_t::eval_type::qpts;
137 const int NC = S_data_t::eval_type::vdim;
138 MFEM_STATIC_ASSERT(qpts == M,
"incompatible dimensions");
139 MFEM_FLOPS_ADD(M*NC);
140 for (
int i = 0; i < M; i++)
142 for (
int j = 0; j < NC; j++)
144 R.val_qpts(i,j,k) *= A[i];
154 template <
int SDim,
int Dim,
typename complex_t>
158 template <
typename complex_t>
164 static const bool uses_Jacobians =
true;
167 static const bool in_values =
false;
168 static const bool in_gradients =
true;
169 static const bool out_values =
false;
170 static const bool out_gradients =
true;
183 template <
typename IR,
typename coeff_t,
int NE>
184 struct CoefficientEval
196 template <
typename T_result_t,
typename Q_t,
typename q_t,
198 static inline MFEM_ALWAYS_INLINE
199 void Action(
const int k,
const T_result_t &F,
200 const Q_t &Q,
const q_t &q, S_data_t &R)
202 const int M = S_data_t::eval_type::qpts;
203 const int NC = S_data_t::eval_type::vdim;
204 MFEM_STATIC_ASSERT(T_result_t::Jt_type::layout_type::dim_1 == M,
205 "incompatible dimensions");
206 MFEM_FLOPS_ADD(M*(1+NC));
207 for (
int i = 0; i < M; i++)
209 const complex_t wi = Q.get(q,i,k) / F.Jt(i,0,0,k);
210 for (
int j = 0; j < NC; j++)
212 R.grad_qpts(i,0,j,k) *= wi;
228 template <
typename T_result_t,
typename Q_t,
typename q_t,
typename asm_type>
229 static inline MFEM_ALWAYS_INLINE
231 const Q_t &Q,
const q_t &q, asm_type &A)
233 const int M = T_result_t::Jt_type::layout_type::dim_1;
234 MFEM_STATIC_ASSERT(asm_type::layout_type::dim_1 == M,
235 "incompatible dimensions");
237 for (
int i = 0; i < M; i++)
240 A[i] = Q.get(q,i,k) / F.Jt(i,0,0,k);
250 template <
int qpts,
typename S_data_t>
251 static inline MFEM_ALWAYS_INLINE
255 const int M = S_data_t::eval_type::qpts;
256 const int NC = S_data_t::eval_type::vdim;
257 MFEM_STATIC_ASSERT(qpts == M,
"incompatible dimensions");
258 MFEM_FLOPS_ADD(M*NC);
259 for (
int i = 0; i < M; i++)
261 for (
int j = 0; j < NC; j++)
263 R.grad_qpts(i,0,j,k) *= A(i,0);
270 template <
typename complex_t>
276 static const bool uses_Jacobians =
true;
279 static const bool in_values =
false;
280 static const bool in_gradients =
true;
281 static const bool out_values =
false;
282 static const bool out_gradients =
true;
296 template <
typename IR,
typename coeff_t,
int NE>
297 struct CoefficientEval
309 template <
typename T_result_t,
typename Q_t,
typename q_t,
311 static inline MFEM_ALWAYS_INLINE
312 void Action(
const int k,
const T_result_t &F,
313 const Q_t &Q,
const q_t &q, S_data_t &R)
315 const int M = S_data_t::eval_type::qpts;
316 const int NC = S_data_t::eval_type::vdim;
317 MFEM_STATIC_ASSERT(T_result_t::Jt_type::layout_type::dim_1 == M,
318 "incompatible dimensions");
319 MFEM_FLOPS_ADD(M*(4+NC*14));
320 for (
int i = 0; i < M; i++)
322 typedef typename T_result_t::Jt_type::data_type
real_t;
323 const real_t J11 = F.Jt(i,0,0,k);
324 const real_t J12 = F.Jt(i,1,0,k);
325 const real_t J21 = F.Jt(i,0,1,k);
326 const real_t J22 = F.Jt(i,1,1,k);
327 const complex_t w_det_J = Q.get(q,i,k) / (J11 * J22 - J21 * J12);
328 for (
int j = 0; j < NC; j++)
330 const complex_t x1 = R.grad_qpts(i,0,j,k);
331 const complex_t x2 = R.grad_qpts(i,1,j,k);
333 const complex_t z1 = J22 * x1 - J21 * x2;
334 const complex_t z2 = J11 * x2 - J12 * x1;
335 R.grad_qpts(i,0,j,k) = w_det_J * (J22 * z1 - J12 * z2);
336 R.grad_qpts(i,1,j,k) = w_det_J * (J11 * z2 - J21 * z1);
352 template <
typename T_result_t,
typename Q_t,
typename q_t,
typename asm_type>
353 static inline MFEM_ALWAYS_INLINE
355 const Q_t &Q,
const q_t &q, asm_type &A)
357 typedef typename T_result_t::Jt_type::data_type
real_t;
358 const int M = T_result_t::Jt_type::layout_type::dim_1;
359 MFEM_STATIC_ASSERT(asm_type::layout_type::dim_1 == M,
360 "incompatible dimensions");
361 MFEM_FLOPS_ADD(16*M);
362 const bool Symm = (asm_type::layout_type::rank == 2);
363 for (
int i = 0; i < M; i++)
365 const real_t J11 = F.Jt(i,0,0,k);
366 const real_t J12 = F.Jt(i,1,0,k);
367 const real_t J21 = F.Jt(i,0,1,k);
368 const real_t J22 = F.Jt(i,1,1,k);
369 const complex_t w_det_J = Q.get(q,i,k) / (J11 * J22 - J21 * J12);
370 internal::MatrixOps<2,2>::Symm<Symm>::Set(
372 + w_det_J * (J12*J12 + J22*J22),
373 - w_det_J * (J11*J12 + J21*J22),
374 + w_det_J * (J11*J11 + J21*J21)
385 template <
int qpts,
typename S_data_t>
386 static inline MFEM_ALWAYS_INLINE
390 const int M = S_data_t::eval_type::qpts;
391 const int NC = S_data_t::eval_type::vdim;
392 MFEM_STATIC_ASSERT(qpts == M,
"incompatible dimensions");
393 MFEM_FLOPS_ADD(6*M*NC);
394 for (
int i = 0; i < M; i++)
396 const complex_t A11 = A(i,0);
397 const complex_t A21 = A(i,1);
398 const complex_t A22 = A(i,2);
399 for (
int j = 0; j < NC; j++)
401 const complex_t x1 = R.grad_qpts(i,0,j,k);
402 const complex_t x2 = R.grad_qpts(i,1,j,k);
403 R.grad_qpts(i,0,j,k) = A11 * x1 + A21 * x2;
404 R.grad_qpts(i,1,j,k) = A21 * x1 + A22 * x2;
411 template <
typename complex_t>
417 static const bool uses_Jacobians =
true;
420 static const bool in_values =
false;
421 static const bool in_gradients =
true;
422 static const bool out_values =
false;
423 static const bool out_gradients =
true;
437 template <
typename IR,
typename coeff_t,
int NE>
438 struct CoefficientEval
450 template <
typename T_result_t,
typename Q_t,
typename q_t,
452 static inline MFEM_ALWAYS_INLINE
453 void Action(
const int k,
const T_result_t &F,
454 const Q_t &Q,
const q_t &q, S_data_t &R)
456 const int M = S_data_t::eval_type::qpts;
457 const int NC = S_data_t::eval_type::vdim;
458 MFEM_STATIC_ASSERT(T_result_t::Jt_type::layout_type::dim_1 == M,
459 "incompatible dimensions");
461 for (
int i = 0; i < M; i++)
463 typedef typename T_result_t::Jt_type::data_type
real_t;
465 const complex_t w_det_J =
467 TAdjDet<real_t>(F.Jt.layout.ind14(i,k).transpose_12(), F.Jt,
470 sMult_AB<false>(adj_J.
layout.transpose_12(), adj_J,
471 R.grad_qpts.
layout.ind14(i,k), R.grad_qpts,
474 sMult_AB<false>(adj_J.
layout, adj_J,
476 R.grad_qpts.
layout.ind14(i,k), R.grad_qpts);
491 template <
typename T_result_t,
typename Q_t,
typename q_t,
typename asm_type>
492 static inline MFEM_ALWAYS_INLINE
494 const Q_t &Q,
const q_t &q, asm_type &A)
496 typedef typename T_result_t::Jt_type::data_type
real_t;
497 const int M = T_result_t::Jt_type::layout_type::dim_1;
498 MFEM_STATIC_ASSERT(asm_type::layout_type::dim_1 == M,
499 "incompatible dimensions");
500 MFEM_FLOPS_ADD(37*M);
501 const bool Symm = (asm_type::layout_type::rank == 2);
502 for (
int i = 0; i < M; i++)
507 TAdjDet<real_t>(F.Jt.layout.ind14(i,k).transpose_12(), F.Jt,
509 internal::MatrixOps<3,3>::Symm<Symm>::Set(
511 u*(B(0,0)*B(0,0)+B(0,1)*B(0,1)+B(0,2)*B(0,2)),
512 u*(B(0,0)*B(1,0)+B(0,1)*B(1,1)+B(0,2)*B(1,2)),
513 u*(B(0,0)*B(2,0)+B(0,1)*B(2,1)+B(0,2)*B(2,2)),
514 u*(B(1,0)*B(1,0)+B(1,1)*B(1,1)+B(1,2)*B(1,2)),
515 u*(B(1,0)*B(2,0)+B(1,1)*B(2,1)+B(1,2)*B(2,2)),
516 u*(B(2,0)*B(2,0)+B(2,1)*B(2,1)+B(2,2)*B(2,2))
527 template <
int qpts,
typename S_data_t>
528 static inline MFEM_ALWAYS_INLINE
532 const int M = S_data_t::eval_type::qpts;
533 const int NC = S_data_t::eval_type::vdim;
534 MFEM_STATIC_ASSERT(qpts == M,
"incompatible dimensions");
535 MFEM_FLOPS_ADD(15*M*NC);
536 for (
int i = 0; i < M; i++)
538 const complex_t A11 = A(i,0);
539 const complex_t A21 = A(i,1);
540 const complex_t A31 = A(i,2);
541 const complex_t A22 = A(i,3);
542 const complex_t A32 = A(i,4);
543 const complex_t A33 = A(i,5);
544 for (
int j = 0; j < NC; j++)
546 const complex_t x1 = R.grad_qpts(i,0,j,k);
547 const complex_t x2 = R.grad_qpts(i,1,j,k);
548 const complex_t x3 = R.grad_qpts(i,2,j,k);
549 R.grad_qpts(i,0,j,k) = A11*x1 + A21*x2 + A31*x3;
550 R.grad_qpts(i,1,j,k) = A21*x1 + A22*x2 + A32*x3;
551 R.grad_qpts(i,2,j,k) = A31*x1 + A32*x2 + A33*x3;
559 #endif // MFEM_TEMPLATE_BILININTEG
static const bool uses_Jacobians
static const bool in_values
static const layout_type layout
static MFEM_ALWAYS_INLINE void Action(const int k, const T_result_t &F, const Q_t &Q, const q_t &q, S_data_t &R)
IntRuleCoefficient< IR, coeff_t, NE >::Type Type
TTensor3< qpts, 2, 2, complex_t > type
TTensor3< qpts, 1, 1, complex_t > type
static const bool out_gradients
TMatrix< qpts, 3, complex_t > type
static MFEM_ALWAYS_INLINE void Assemble(const int k, const T_result_t &F, const Q_t &Q, const q_t &q, asm_type &A)
static MFEM_ALWAYS_INLINE void Assemble(const int k, const T_result_t &F, const Q_t &Q, const q_t &q, asm_type &A)
static MFEM_ALWAYS_INLINE void MultAssembled(const int k, const TMatrix< qpts, 3, complex_t > &A, S_data_t &R)
IntRuleCoefficient< IR, coeff_t, NE >::Type Type
kernel_t< SDim, Dim, complex_t > type
static const bool in_gradients
TVector< qpts, complex_t > type
static MFEM_ALWAYS_INLINE void Action(const int k, const T_result_t &F, const Q_t &Q, const q_t &q, S_data_t &R)
TMatrix< qpts, 6, complex_t > type
static MFEM_ALWAYS_INLINE void Assemble(const int k, const T_result_t &F, const Q_t &Q, const q_t &q, asm_type &A)
static MFEM_ALWAYS_INLINE void MultAssembled(const int k, const TMatrix< qpts, 1, complex_t > &A, S_data_t &R)
static MFEM_ALWAYS_INLINE void Assemble(const int k, const T_result_t &F, const Q_t &Q, const q_t &q, TVector< qpts, complex_t > &A)
IntRuleCoefficient< IR, coeff_t, NE >::Type Type
IntRuleCoefficient< IR, coeff_t, NE >::Type Type
static MFEM_ALWAYS_INLINE void Action(const int k, const T_result_t &F, const Q_t &Q, const q_t &q, S_data_t &R)
TMatrix< qpts, 1, complex_t > type
static MFEM_ALWAYS_INLINE void MultAssembled(const int k, const TMatrix< qpts, 6, complex_t > &A, S_data_t &R)
TVector< qpts, complex_t > type
static const bool out_values
static MFEM_ALWAYS_INLINE void Action(const int k, const T_result_t &F, const Q_t &Q, const q_t &q, S_data_t &R)
static MFEM_ALWAYS_INLINE void MultAssembled(const int k, const TVector< qpts, complex_t > &A, S_data_t &R)
TTensor3< qpts, 3, 3, complex_t > type
void Scale(const data_t scale)
TIntegrator(const coefficient_type &c)