21#include "rheolef/fem_on_pointset.h"
22#include "rheolef/piola_util.h"
47template <
class T,
class M,
class Value, details::differentiate_option::type Diff>
48struct evaluate_internal {
53 const details::differentiate_option& gopt,
54 Eigen::Matrix<Value,Eigen::Dynamic,Eigen::Dynamic>& value)
59template <
class T,
class M,
class Value, details::differentiate_option::type Diff>
60struct evaluate_on_side_internal {
62 const fem_on_pointset_rep<T>& obj,
63 const geo_basic<T,M>& omega_K,
65 const side_information_type& sid,
66 const details::differentiate_option& gopt,
67 Eigen::Matrix<Value,Eigen::Dynamic,Eigen::Dynamic>& value)
77template<
class M,
class Value>
82 const Eigen::Matrix<Value,Eigen::Dynamic,Eigen::Dynamic>& hat_phij_xi,
83 Eigen::Matrix<Value,Eigen::Dynamic,Eigen::Dynamic>& value)
const
85 const piola_fem<T>& pf = _bops.get_basis().get_piola_fem();
86 if (! pf.transform_need_piola()) {
92 const Eigen::Matrix<piola<T>,Eigen::Dynamic,1>&
piola = _pops.get_piola (omega_K, K);
95 value.resize (loc_nnod, loc_ndof);
96 for (
size_type loc_inod = 0; loc_inod < loc_nnod; ++loc_inod) {
97 for (
size_type loc_jdof = 0; loc_jdof < loc_ndof; ++loc_jdof) {
98 const Value& hat_u = hat_phij_xi (loc_inod, loc_jdof);
99 Value&
u = value (loc_inod, loc_jdof);
100 pf.transform (
piola[loc_inod], hat_u,
u);
104template <
class T,
class M,
class Value>
105struct evaluate_internal<
T,
M,Value,details::differentiate_option::none> {
110 const details::differentiate_option& gopt,
111 Eigen::Matrix<Value,Eigen::Dynamic,Eigen::Dynamic>& value)
const
113 const Eigen::Matrix<Value,Eigen::Dynamic,Eigen::Dynamic>& hat_phij_xi = obj.
_bops.template evaluate<Value> (K);
117template <
class T,
class M,
class Value>
118struct evaluate_on_side_internal<
T,
M,Value,details::differentiate_option::none> {
120 const fem_on_pointset_rep<T>& obj,
121 const geo_basic<T,M>& omega_K,
122 const geo_element& K,
123 const side_information_type& sid,
124 const details::differentiate_option& gopt,
125 Eigen::Matrix<Value,Eigen::Dynamic,Eigen::Dynamic>& value)
const
127 const Eigen::Matrix<Value,Eigen::Dynamic,Eigen::Dynamic>& hat_phij_xi = obj._bops.template evaluate_on_side<Value> (K, sid);
128 obj._evaluate_post_piola (omega_K, K, hat_phij_xi, value);
136template<
class T,
class M,
class Value,
class GradValue>
139grad_evaluate_post_piola (
140 const fem_on_pointset_rep<T>& obj,
141 const geo_basic<T,M>& omega_K,
142 const geo_element& K,
143 const details::differentiate_option& gopt,
144 const Eigen::Matrix<Value,Eigen::Dynamic,Eigen::Dynamic>&
146 const Eigen::Matrix<GradValue,Eigen::Dynamic,Eigen::Dynamic>&
148 Eigen::Matrix<GradValue,Eigen::Dynamic,Eigen::Dynamic>& value)
151 reference_element hat_K = K;
152 const Eigen::Matrix<piola<T>,Eigen::Dynamic,1>& piola = obj._pops.get_piola (omega_K, K);
153 const piola_fem<T>& pf = obj._bops.get_basis().get_piola_fem();
154 size_type loc_nnod = hat_grad_phij_xi.rows();
155 size_type loc_ndof = hat_grad_phij_xi.cols();
156 value.resize (loc_nnod,loc_ndof);
157 for (
size_type loc_inod = 0; loc_inod < loc_nnod; ++loc_inod) {
158 for (
size_type loc_jdof = 0; loc_jdof < loc_ndof; ++loc_jdof) {
159 const Value& hat_u = hat_phij_xi (loc_inod, loc_jdof);
160 const GradValue& hat_grad_u = hat_grad_phij_xi (loc_inod, loc_jdof);
161 GradValue&
grad_u = value (loc_inod, loc_jdof);
162 pf.grad_transform (piola[loc_inod], hat_u, hat_grad_u, gopt,
grad_u);
166template<
class T,
class M>
169grad_evaluate_post_piola (
170 const fem_on_pointset_rep<T>& obj,
171 const geo_basic<T,M>& omega_K,
172 const geo_element& K,
173 const details::differentiate_option& gopt,
174 const Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic>&
176 Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic>& value)
178 fatal_macro(
"grad_evaluate: undefined scalar-valued grad");
180template <
class T,
class M,
class GradValue>
181struct evaluate_internal<
T,
M,GradValue,details::differentiate_option::gradient> {
183 const fem_on_pointset_rep<T>& obj,
184 const geo_basic<T,M>& omega_K,
185 const geo_element& K,
186 const details::differentiate_option& gopt,
187 Eigen::Matrix<GradValue,Eigen::Dynamic,Eigen::Dynamic>& value)
189 reference_element hat_K = K;
191 const Eigen::Matrix<Value,Eigen::Dynamic,Eigen::Dynamic>&
192 hat_phij_xi = obj._bops.template evaluate<Value> (hat_K);
193 const Eigen::Matrix<GradValue,Eigen::Dynamic,Eigen::Dynamic>&
194 hat_grad_phij_xi = obj._bops.template grad_evaluate<GradValue> (hat_K);
195 grad_evaluate_post_piola (obj, omega_K, K, gopt, hat_phij_xi, hat_grad_phij_xi, value);
198template <
class T,
class M>
199struct evaluate_internal<
T,
M,
T,details::differentiate_option::gradient> {
201 const fem_on_pointset_rep<T>& obj,
202 const geo_basic<T,M>& omega_K,
203 const geo_element& K,
204 const details::differentiate_option& gopt,
205 Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic>& value)
210template <
class T,
class M,
class GradValue>
211struct evaluate_on_side_internal<
T,
M,GradValue,details::differentiate_option::gradient> {
213 const fem_on_pointset_rep<T>& obj,
214 const geo_basic<T,M>& omega_K,
215 const geo_element& K,
216 const side_information_type& sid,
217 const details::differentiate_option& gopt,
218 Eigen::Matrix<GradValue,Eigen::Dynamic,Eigen::Dynamic>& value)
220 reference_element hat_K = K;
222 const Eigen::Matrix<Value,Eigen::Dynamic,Eigen::Dynamic>&
223 hat_phij_xi = obj._bops.template evaluate_on_side<Value> (hat_K,sid);
224 const Eigen::Matrix<GradValue,Eigen::Dynamic,Eigen::Dynamic>&
225 hat_grad_phij_xi = obj._bops.template grad_evaluate_on_side<GradValue> (hat_K, sid);
226 grad_evaluate_post_piola (obj, omega_K, K, gopt, hat_phij_xi, hat_grad_phij_xi, value);
229template <
class T,
class M>
230struct evaluate_on_side_internal<
T,
M,
T,details::differentiate_option::gradient> {
232 const fem_on_pointset_rep<T>& obj,
233 const geo_basic<T,M>& omega_K,
234 const geo_element& K,
235 const side_information_type& sid,
236 const details::differentiate_option& gopt,
237 Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic>& value)
239 fatal_macro(
"grad_evaluate_on_side: undefined scalar-valued grad");
247template<
class T,
class M,
class Value,
class GradValue>
250div_evaluate_internal (
251 const fem_on_pointset_rep<T>& obj,
252 const geo_basic<T,M>& omega_K,
253 const geo_element& K,
254 const details::differentiate_option& gopt,
255 const Eigen::Matrix<GradValue,Eigen::Dynamic,Eigen::Dynamic>& grad_value,
256 Eigen::Matrix<Value,Eigen::Dynamic,Eigen::Dynamic>& value)
260 <<
"-valued div operator");
263template <
class T,
class M>
265div_evaluate_internal (
269 const details::differentiate_option& gopt,
270 const Eigen::Matrix<
tensor_basic<T>,Eigen::Dynamic,Eigen::Dynamic>& grad_value,
271 Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic>& value)
274 value.resize (grad_value.rows(), grad_value.cols());
275 for (
size_type loc_inod = 0, loc_nnod = value.rows(); loc_inod < loc_nnod; ++loc_inod) {
276 for (
size_type loc_jdof = 0, loc_ndof = value.cols(); loc_jdof < loc_ndof; ++loc_jdof) {
277 value (loc_inod,loc_jdof) =
tr (grad_value(loc_inod,loc_jdof));
280template <
class T,
class M,
class Value>
281struct evaluate_internal<
T,
M,Value,details::differentiate_option::divergence> {
283 const fem_on_pointset_rep<T>& obj,
284 const geo_basic<T,M>& omega_K,
285 const geo_element& K,
286 const details::differentiate_option& gopt,
287 Eigen::Matrix<Value,Eigen::Dynamic,Eigen::Dynamic>& value)
const
291 Eigen::Matrix<GradValue,Eigen::Dynamic,Eigen::Dynamic> grad_value;
292 obj.template evaluate<M,GradValue,details::differentiate_option::gradient> (omega_K, K, gopt, grad_value);
293 div_evaluate_internal (obj, omega_K, K, gopt, grad_value, value);
296template <
class T,
class M,
class Value>
297struct evaluate_on_side_internal<
T,
M,Value,details::differentiate_option::divergence> {
299 const fem_on_pointset_rep<T>& obj,
300 const geo_basic<T,M>& omega_K,
301 const geo_element& K,
302 const side_information_type& sid,
303 const details::differentiate_option& gopt,
304 Eigen::Matrix<Value,Eigen::Dynamic,Eigen::Dynamic>& value)
const
308 Eigen::Matrix<GradValue,Eigen::Dynamic,Eigen::Dynamic> grad_value;
309 obj.template evaluate_on_side<M,GradValue,details::differentiate_option::gradient> (omega_K, K, sid, gopt, grad_value);
310 div_evaluate_internal (obj, omega_K, K, gopt, grad_value, value);
318template<
class T,
class M,
class Value>
321curl_evaluate_internal (
322 const fem_on_pointset_rep<T>& obj,
323 const geo_basic<T,M>& omega_K,
324 const geo_element& K,
325 const details::differentiate_option& gopt,
326 Eigen::Matrix<Value,Eigen::Dynamic,Eigen::Dynamic>& value)
329 <<
"-valued curl operator");
332template <
class T,
class M>
334curl_evaluate_internal (
338 const details::differentiate_option& gopt,
339 Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic>& value)
342 Eigen::Matrix<tensor_basic<T>,Eigen::Dynamic,Eigen::Dynamic> value_grad;
344 value.resize (value_grad.rows(), value_grad.cols());
345 for (
size_type loc_inod = 0, loc_nnod = value.rows(); loc_inod < loc_nnod; ++loc_inod) {
346 for (
size_type loc_jdof = 0, loc_ndof = value.cols(); loc_jdof < loc_ndof; ++loc_jdof) {
348 value (loc_inod,loc_jdof) =
g(1,0) -
g(0,1);
354template <
class T,
class M>
356curl_evaluate_internal (
360 const details::differentiate_option& gopt,
361 Eigen::Matrix<
point_basic<T>,Eigen::Dynamic,Eigen::Dynamic>& value)
365 check_macro (
d >= 2,
"curl: unexpected dimension d="<<
d<<
" for mesh \""<<omega_K.name()<<
"\"");
369 Eigen::Matrix<point_basic<T>,Eigen::Dynamic,Eigen::Dynamic> value_grad;
371 value.resize (value_grad.rows(), value_grad.cols());
372 for (
size_type loc_inod = 0, loc_nnod = value.rows(); loc_inod < loc_nnod; ++loc_inod) {
373 for (
size_type loc_jdof = 0, loc_ndof = value.cols(); loc_jdof < loc_ndof; ++loc_jdof) {
389 const Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic>&
390 hat_phij_xi = obj.
_bops.template evaluate<T> (hat_K);
391 const Eigen::Matrix<piola<T>,Eigen::Dynamic,1>& piola = obj.
_pops.get_piola (omega_K, K);
394 for (
size_type loc_inod = 0; loc_inod < loc_nnod; ++loc_inod) {
395 T r = piola[loc_inod].F[i_comp_r];
396 check_macro (1+abs(r) != 1,
"curl(): singular axisymmetric (1/r) weight (HINT: avoid interpolate() or change quadrature formulae)");
398 for (
size_type loc_jdof = 0; loc_jdof < loc_ndof; ++loc_jdof) {
399 const T& w = hat_phij_xi (loc_inod,loc_jdof);
410 Eigen::Matrix<tensor_basic<T>,Eigen::Dynamic,Eigen::Dynamic> value_grad;
412 value.resize (value_grad.rows(), value_grad.cols());
413 for (
size_type loc_inod = 0, loc_nnod = value.rows(); loc_inod < loc_nnod; ++loc_inod) {
414 for (
size_type loc_jdof = 0, loc_ndof = value.cols(); loc_jdof < loc_ndof; ++loc_jdof) {
417 c[0] =
g(2,1) -
g(1,2);
418 c[1] =
g(0,2) -
g(2,0);
419 c[2] =
g(1,0) -
g(0,1);
423template <
class T,
class M,
class Value>
424struct evaluate_internal<
T,
M,Value,details::differentiate_option::curl> {
426 const fem_on_pointset_rep<T>& obj,
427 const geo_basic<T,M>& omega_K,
428 const geo_element& K,
429 const details::differentiate_option& gopt,
430 Eigen::Matrix<Value,Eigen::Dynamic,Eigen::Dynamic>& value)
const
432 curl_evaluate_internal (obj, omega_K, K, gopt, value);
435template <
class T,
class M,
class Value>
436struct evaluate_on_side_internal<
T,
M,Value,details::differentiate_option::curl> {
438 const fem_on_pointset_rep<T>& obj,
439 const geo_basic<T,M>& omega_K,
440 const geo_element& K,
441 const side_information_type& sid,
442 const details::differentiate_option& gopt,
443 Eigen::Matrix<Value,Eigen::Dynamic,Eigen::Dynamic>& value)
const
452template <
class M,
class Value, details::differentiate_option::type Diff>
458 Eigen::Matrix<Value,Eigen::Dynamic,Eigen::Dynamic>& value)
const
460 evaluate_internal<T,M,Value,Diff> eval;
461 eval (*
this, omega_K, K, gopt, value);
464template <
class M,
class Value, details::differentiate_option::type Diff>
471 Eigen::Matrix<Value,Eigen::Dynamic,Eigen::Dynamic>& value)
const
473 evaluate_on_side_internal<T,M,Value,Diff> eval;
474 eval (*
this, omega_K, K, sid, gopt, value);
479#define _RHEOLEF_instanciation(T) \
480template class fem_on_pointset_rep<T>; \
482#define _RHEOLEF_instanciation_value(T,M,Value,Diff) \
483template void fem_on_pointset_rep<T>::evaluate<M,Value,Diff> ( \
484 const geo_basic<T,M>& omega_K, \
485 const geo_element& K, \
486 const details::differentiate_option& gopt, \
487 Eigen::Matrix<Value,Eigen::Dynamic,Eigen::Dynamic>& value) const; \
488template void fem_on_pointset_rep<T>::evaluate_on_side<M,Value,Diff> ( \
489 const geo_basic<T,M>& omega_K, \
490 const geo_element& K, \
491 const side_information_type& sid, \
492 const details::differentiate_option& gopt, \
493 Eigen::Matrix<Value,Eigen::Dynamic,Eigen::Dynamic>& value) const; \
495#define _RHEOLEF_instanciation_values(T,M,Diff) \
496_RHEOLEF_instanciation_value(T,M,T,Diff) \
497_RHEOLEF_instanciation_value(T,M,point_basic<T>,Diff) \
498_RHEOLEF_instanciation_value(T,M,tensor_basic<T>,Diff) \
499_RHEOLEF_instanciation_value(T,M,tensor3_basic<T>,Diff) \
500_RHEOLEF_instanciation_value(T,M,tensor4_basic<T>,Diff) \
502#define _RHEOLEF_instanciation_evaluate(T,M) \
503_RHEOLEF_instanciation_values(T,M,details::differentiate_option::none) \
504_RHEOLEF_instanciation_values(T,M,details::differentiate_option::gradient) \
505_RHEOLEF_instanciation_values(T,M,details::differentiate_option::divergence) \
506_RHEOLEF_instanciation_values(T,M,details::differentiate_option::curl) \
510#ifdef _RHEOLEF_HAVE_MPI
#define _RHEOLEF_instanciation(T, M, A)
field::size_type size_type
see the Float page for the full documentation
basis_on_pointset< T > _bops
void _evaluate_post_piola(const geo_basic< T, M > &omega_K, const geo_element &K, const Eigen::Matrix< Value, Eigen::Dynamic, Eigen::Dynamic > &hat_phij_xi, Eigen::Matrix< Value, Eigen::Dynamic, Eigen::Dynamic > &value) const
reference_element::size_type size_type
void evaluate(const geo_basic< T, M > &omega_K, const geo_element &K, const details::differentiate_option &gopt, Eigen::Matrix< Value, Eigen::Dynamic, Eigen::Dynamic > &value) const
void initialize(const basis_basic< T > &fem_basis, const piola_on_pointset< T > &pops)
void evaluate_on_side(const geo_basic< T, M > &omega_K, const geo_element &K, const side_information_type &sid, const details::differentiate_option &gopt, Eigen::Matrix< Value, Eigen::Dynamic, Eigen::Dynamic > &value) const
piola_on_pointset< T > _pops
generic mesh with rerefence counting
see the geo_element page for the full documentation
bool has_quadrature() const
const quadrature< T > & get_quadrature() const
const basis_basic< T > & get_nodal_basis() const
see the reference_element page for the full documentation
#define fatal_macro(message)
#define _RHEOLEF_instanciation_evaluate(T, M)
check_macro(expr1.have_homogeneous_space(Xh1), "dual(expr1,expr2); expr1 should have homogeneous space. HINT: use dual(interpolate(Xh, expr1),expr2)")
const std::string & valued_name(valued_type valued_tag)
This file is part of Rheolef.
U tr(const tensor_basic< U > &a, size_t d=3)
undeterminated_basic< typename scalar_traits< T >::type > type
point_basic< typename scalar_traits< T >::type > type
static const valued_type value