Medusa  1.1
Coordinate Free Mehless Method implementation
ExplicitOperators.hpp
Go to the documentation of this file.
1 #ifndef MEDUSA_BITS_OPERATORS_EXPLICITOPERATORS_HPP_
2 #define MEDUSA_BITS_OPERATORS_EXPLICITOPERATORS_HPP_
3 
12 #include "ShapeStorage_fwd.hpp"
13 
14 namespace mm {
15 
21 #define NODE_ASSERTS(node) \
22  assert_msg(0 <= (node) && (node) < ss->size(), "Node index %d must be in range " \
23  "[0, num_nodes = %d).", ss->size());
24 
26 template <typename shape_storage_type>
27 template <typename op_family_t, typename scalar_field_t>
29 ExplicitOperators<shape_storage_type>::apply(const scalar_field_t& u, int node,
30  typename op_family_t::operator_t o) const {
31  NODE_ASSERTS(node);
32  int op_idx = op_family_t::index(o);
33  typename scalar_type<scalar_field_t>::type ret = 0.0;
34  for (int j = 0; j < ss->supportSize(node); ++j)
35  ret += ss->template get<op_family_t>(op_idx, node, j) * u[ss->support(node, j)];
36  return ret;
37 }
38 
39 template <typename shape_storage_type>
40 template <typename scalar_field_t>
41 Vec<typename scalar_type<scalar_field_t>::type, ExplicitOperators<shape_storage_type>::dim>
42 ExplicitOperators<shape_storage_type>::grad(const scalar_field_t& u, int node) const {
43  NODE_ASSERTS(node);
44  Vec<typename scalar_type<scalar_field_t>::type, ExplicitOperators<shape_storage_type>::dim> ret;
45  ret.setZero();
46  for (int var = 0; var < dim; ++var) {
47  for (int i = 0; i < ss->supportSize(node); ++i) {
48  ret[var] += ss->d1(var, node, i) * u[ss->support(node, i)];
49  }
50  }
51  return ret;
52 }
53 
54 
55 template <typename shape_storage_type>
56 template <typename scalar_field_t>
59  const scalar_field_t& u, int node, const vector_t& normal,
60  typename scalar_type<scalar_field_t>::type val) const {
61  NODE_ASSERTS(node);
62  assert_msg(ss->support(node, 0) == node, "First support node should be the node itself.");
63  scalar_t denominator = 0;
64  for (int d = 0; d < dim; ++d) {
65  for (int i = 1; i < ss->supportSize(node); ++i) {
66  val -= normal[d] * ss->d1(d, node, i) * u[ss->support(node, i)];
67  }
68  // i = 0
69  denominator += normal[d] * ss->d1(d, node, 0);
70  }
71  assert_msg(std::abs(denominator) >= 1e-14,
72  "Node %d has no effect on the flux in direction %s. The cause of this might be wrong"
73  " normal direction, bad neighbourhood choice or bad nodal positions.", node, normal);
74  return val / denominator;
75 }
77 
79 template <typename S>
80 std::ostream& operator<<(std::ostream& os, const ExplicitOperators<S>& op) {
81  if (!op.hasShapes()) {
82  return os << "Explicit operators without any linked storage.";
83  }
84  return os << "Explicit operators over storage: " << *op.ss;
85 }
86 
87 template <typename Derived, typename vec_t, typename OpFamilies>
88 ExplicitOperators<Derived>
90  return ExplicitOperators<Derived>(*static_cast<const Derived*>(this));
91 }
92 
93 } // namespace mm
94 
95 #endif // MEDUSA_BITS_OPERATORS_EXPLICITOPERATORS_HPP_
mm::ExplicitOperators::ss
const shape_storage_t * ss
Pointer to shape storage, but name is shortened for readability.
Definition: ExplicitOperators_fwd.hpp:45
mm
Root namespace for the whole library.
Definition: Gaussian.hpp:14
scalar_t
Scalar scalar_t
Type of the elements, alias of Scalar.
Definition: MatrixBaseAddons.hpp:16
mm::ExplicitOperators::neumann
scalar_type< scalar_field_t >::type neumann(const scalar_field_t &u, int node, const vector_t &normal, typename scalar_type< scalar_field_t >::type val) const
Calculates a new field value a for u at node, such that the Neumann boundary condition du/dn(node) = ...
ShapeStorage_fwd.hpp
dim
@ dim
Number of elements of this matrix.
Definition: MatrixBaseAddons.hpp:14
mm::operator<<
std::ostream & operator<<(std::ostream &os, const Gaussian< S > &b)
Output basic information about given Gaussian RBF.
Definition: Gaussian.hpp:37
assert_msg
#define assert_msg(cond,...)
Assert with better error reporting.
Definition: assert.hpp:75
NODE_ASSERTS
#define NODE_ASSERTS(node)
Performs all asserts for a given node: asserts that node index is valid and that shape functions were...
Definition: ExplicitOperators.hpp:21
mm::ExplicitOperators
A class for evaluating typical operators needed in spatial discretization.
Definition: ExplicitOperators_fwd.hpp:35
mm::ExplicitOperators::apply
scalar_type< scalar_field_t >::type apply(const scalar_field_t &u, int node, typename op_family_t::operator_t o) const
Returns an approximation of applying the requested operator to field u at node node.
assert.hpp
mm::ExplicitOperators::grad
Vec< typename scalar_type< scalar_field_t >::type, ExplicitOperators< shape_storage_type >::dim > grad(const scalar_field_t &u, int node) const
Returns an approximation of gradient of field u in node-th node.
Vec.hpp
ExplicitOperators_fwd.hpp
mm::ExplicitOperators::hasShapes
bool hasShapes() const
Returns true if operators have a non-null pointer to storage and false otherwise.
Definition: ExplicitOperators_fwd.hpp:69
mm::scalar_type::type
scalar_field_t::Scalar type
Default scalar type.
Definition: traits.hpp:19
mm::ExplicitOperators::dim
@ dim
Dimensionality of the function domain.
Definition: ExplicitOperators_fwd.hpp:41
mm::ShapeStorage::explicitOperators
ExplicitOperators< Derived > explicitOperators() const
Construct explicit operators over this storage.
Definition: ExplicitOperators.hpp:89