DUNE-FEM (unstable)

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// SPDX-FileCopyrightInfo: Copyright © DUNE Project contributors, see file LICENSE.md in module root
4// SPDX-License-Identifier: LicenseRef-GPL-2.0-only-with-DUNE-exception
5#ifndef DUNE_DUAL_Q1_LOCALBASIS_HH
6#define DUNE_DUAL_Q1_LOCALBASIS_HH
7
8#include <array>
9#include <numeric>
10
13
14#include <dune/localfunctions/common/localbasis.hh>
15
16namespace Dune
17{
28 template<class D, class R, int dim>
30 {
31 public:
34
35 void setCoefficients(const std::array<Dune::FieldVector<R, (1<<dim)> ,(1<<dim)>& coefficients)
36 {
37 coefficients_ = coefficients;
38 }
39
41 unsigned int size () const
42 {
43 return 1<<dim;
44 }
45
47 inline void evaluateFunction (const typename Traits::DomainType& in,
48 std::vector<typename Traits::RangeType>& out) const
49 {
50 // compute q1 values
51 std::vector<typename Traits::RangeType> q1Values(size());
52
53 for (size_t i=0; i<size(); i++) {
54
55 q1Values[i] = 1;
56
57 for (int j=0; j<dim; j++)
58 // if j-th bit of i is set multiply with in[j], else with 1-in[j]
59 q1Values[i] *= (i & (1<<j)) ? in[j] : 1-in[j];
60
61 }
62
63 // compute the dual values by using that they are linear combinations of q1 functions
64 out.resize(size());
65 for (size_t i=0; i<size(); i++)
66 out[i] = 0;
67
68 for (size_t i=0; i<size(); i++)
69 for (size_t j=0; j<size(); j++)
70 out[i] += coefficients_[i][j]*q1Values[j];
71
72
73 }
74
76 inline void
77 evaluateJacobian (const typename Traits::DomainType& in, // position
78 std::vector<typename Traits::JacobianType>& out) const // return value
79 {
80 // compute q1 jacobians
81 std::vector<typename Traits::JacobianType> q1Jacs(size());
82
83 // Loop over all shape functions
84 for (size_t i=0; i<size(); i++) {
85
86 // Loop over all coordinate directions
87 for (int j=0; j<dim; j++) {
88
89 // Initialize: the overall expression is a product
90 // if j-th bit of i is set to -1, else 1
91 q1Jacs[i][0][j] = (i & (1<<j)) ? 1 : -1;
92
93 for (int k=0; k<dim; k++) {
94
95 if (j!=k)
96 // if k-th bit of i is set multiply with in[j], else with 1-in[j]
97 q1Jacs[i][0][j] *= (i & (1<<k)) ? in[k] : 1-in[k];
98
99 }
100
101 }
102
103 }
104
105 // compute the dual jacobians by using that they are linear combinations of q1 functions
106 out.resize(size());
107 for (size_t i=0; i<size(); i++)
108 out[i] = 0;
109
110 for (size_t i=0; i<size(); i++)
111 for (size_t j=0; j<size(); j++)
112 out[i].axpy(coefficients_[i][j],q1Jacs[j]);
113
114 }
115
117 void partial (const std::array<unsigned int, dim>& order,
118 const typename Traits::DomainType& in, // position
119 std::vector<typename Traits::RangeType>& out) const // return value
120 {
121 auto totalOrder = std::accumulate(order.begin(), order.end(), 0);
122 if (totalOrder == 0) {
123 evaluateFunction(in, out);
124 } else {
125 DUNE_THROW(NotImplemented, "Desired derivative order is not implemented");
126 }
127 }
128
130 unsigned int order () const
131 {
132 return 1;
133 }
134
135 private:
136 std::array<Dune::FieldVector<R, (1<<dim)> ,(1<<dim)> coefficients_;
137 };
138}
139#endif
Dual Lagrange shape functions of order 1 on the reference cube.
Definition: dualq1localbasis.hh:30
unsigned int size() const
number of shape functions
Definition: dualq1localbasis.hh:41
unsigned int order() const
Polynomial order of the shape functions.
Definition: dualq1localbasis.hh:130
void evaluateFunction(const typename Traits::DomainType &in, std::vector< typename Traits::RangeType > &out) const
Evaluate all shape functions.
Definition: dualq1localbasis.hh:47
void evaluateJacobian(const typename Traits::DomainType &in, std::vector< typename Traits::JacobianType > &out) const
Evaluate Jacobian of all shape functions.
Definition: dualq1localbasis.hh:77
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:117
A dense n x m matrix.
Definition: fmatrix.hh:117
vector space out of a tensor product of fields.
Definition: fvector.hh:91
Default exception for dummy implementations.
Definition: exceptions.hh:263
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:218
constexpr T accumulate(Range &&range, T value, F &&f)
Accumulate values.
Definition: hybridutilities.hh:279
Dune namespace.
Definition: alignedallocator.hh:13
Type traits for LocalBasisVirtualInterface.
Definition: localbasis.hh:35
D DomainType
domain type
Definition: localbasis.hh:43
Creative Commons License   |  Legal Statements / Impressum  |  Hosted by TU Dresden  |  generated with Hugo v0.111.3 (Nov 12, 23:30, 2024)