Dune Core Modules (2.8.0)

dualq1localbasis.hh
1// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
2// vi: set et ts=4 sw=2 sts=2:
3#ifndef DUNE_DUAL_Q1_LOCALBASIS_HH
4#define DUNE_DUAL_Q1_LOCALBASIS_HH
5
6#include <array>
7#include <numeric>
8
11
12#include <dune/localfunctions/common/localbasis.hh>
13
14namespace Dune
15{
25 template<class D, class R, int dim>
27 {
28 public:
31
32 void setCoefficients(const std::array<Dune::FieldVector<R, (1<<dim)> ,(1<<dim)>& coefficients)
33 {
34 coefficients_ = coefficients;
35 }
36
38 unsigned int size () const
39 {
40 return 1<<dim;
41 }
42
44 inline void evaluateFunction (const typename Traits::DomainType& in,
45 std::vector<typename Traits::RangeType>& out) const
46 {
47 // compute q1 values
48 std::vector<typename Traits::RangeType> q1Values(size());
49
50 for (size_t i=0; i<size(); i++) {
51
52 q1Values[i] = 1;
53
54 for (int j=0; j<dim; j++)
55 // if j-th bit of i is set multiply with in[j], else with 1-in[j]
56 q1Values[i] *= (i & (1<<j)) ? in[j] : 1-in[j];
57
58 }
59
60 // compute the dual values by using that they are linear combinations of q1 functions
61 out.resize(size());
62 for (size_t i=0; i<size(); i++)
63 out[i] = 0;
64
65 for (size_t i=0; i<size(); i++)
66 for (size_t j=0; j<size(); j++)
67 out[i] += coefficients_[i][j]*q1Values[j];
68
69
70 }
71
73 inline void
74 evaluateJacobian (const typename Traits::DomainType& in, // position
75 std::vector<typename Traits::JacobianType>& out) const // return value
76 {
77 // compute q1 jacobians
78 std::vector<typename Traits::JacobianType> q1Jacs(size());
79
80 // Loop over all shape functions
81 for (size_t i=0; i<size(); i++) {
82
83 // Loop over all coordinate directions
84 for (int j=0; j<dim; j++) {
85
86 // Initialize: the overall expression is a product
87 // if j-th bit of i is set to -1, else 1
88 q1Jacs[i][0][j] = (i & (1<<j)) ? 1 : -1;
89
90 for (int k=0; k<dim; k++) {
91
92 if (j!=k)
93 // if k-th bit of i is set multiply with in[j], else with 1-in[j]
94 q1Jacs[i][0][j] *= (i & (1<<k)) ? in[k] : 1-in[k];
95
96 }
97
98 }
99
100 }
101
102 // compute the dual jacobians by using that they are linear combinations of q1 functions
103 out.resize(size());
104 for (size_t i=0; i<size(); i++)
105 out[i] = 0;
106
107 for (size_t i=0; i<size(); i++)
108 for (size_t j=0; j<size(); j++)
109 out[i].axpy(coefficients_[i][j],q1Jacs[j]);
110
111 }
112
114 void partial (const std::array<unsigned int, dim>& order,
115 const typename Traits::DomainType& in, // position
116 std::vector<typename Traits::RangeType>& out) const // return value
117 {
118 auto totalOrder = std::accumulate(order.begin(), order.end(), 0);
119 if (totalOrder == 0) {
120 evaluateFunction(in, out);
121 } else {
122 DUNE_THROW(NotImplemented, "Desired derivative order is not implemented");
123 }
124 }
125
127 unsigned int order () const
128 {
129 return 1;
130 }
131
132 private:
133 std::array<Dune::FieldVector<R, (1<<dim)> ,(1<<dim)> coefficients_;
134 };
135}
136#endif
Dual Lagrange shape functions of order 1 on the reference cube.
Definition: dualq1localbasis.hh:27
unsigned int size() const
number of shape functions
Definition: dualq1localbasis.hh:38
unsigned int order() const
Polynomial order of the shape functions.
Definition: dualq1localbasis.hh:127
void evaluateFunction(const typename Traits::DomainType &in, std::vector< typename Traits::RangeType > &out) const
Evaluate all shape functions.
Definition: dualq1localbasis.hh:44
void evaluateJacobian(const typename Traits::DomainType &in, std::vector< typename Traits::JacobianType > &out) const
Evaluate Jacobian of all shape functions.
Definition: dualq1localbasis.hh:74
void partial(const std::array< unsigned int, dim > &order, const typename Traits::DomainType &in, std::vector< typename Traits::RangeType > &out) const
Evaluate partial derivatives of all shape functions.
Definition: dualq1localbasis.hh:114
A dense n x m matrix.
Definition: fmatrix.hh:69
vector space out of a tensor product of fields.
Definition: fvector.hh:95
Default exception for dummy implementations.
Definition: exceptions.hh:261
Implements a matrix constructed from a given type representing a field and compile-time given number ...
Implements a vector constructed from a given type representing a field and a compile-time given size.
#define DUNE_THROW(E, m)
Definition: exceptions.hh:216
constexpr T accumulate(Range &&range, T value, F &&f)
Accumulate values.
Definition: hybridutilities.hh:289
Dune namespace.
Definition: alignedallocator.hh:11
Type traits for LocalBasisVirtualInterface.
Definition: localbasis.hh:32
D DomainType
domain type
Definition: localbasis.hh:43
Creative Commons License   |  Legal Statements / Impressum  |  Hosted by TU Dresden  |  generated with Hugo v0.111.3 (Dec 22, 23:30, 2024)