OGS
NsAndWeight.h
Go to the documentation of this file.
1
11#pragma once
12
13#include <Eigen/Core>
14
18
20{
21
31template <typename ShapeFunction, typename LowerOrderShapeFunction,
32 typename ShapeMatrix, typename LowerOrderShapeMatrix>
34{
35 static_assert(
36 ShapeFunction::ORDER > LowerOrderShapeFunction::ORDER,
37 "Shape function orders must be different from each other. For same "
38 "orders the class template specialization below should be used.");
39
40 static_assert(ShapeFunction::ORDER == 2);
41 static_assert(LowerOrderShapeFunction::ORDER < 2);
42
43 NsAndWeight(ShapeMatrix&& N1, LowerOrderShapeMatrix&& N2,
44 double const weight)
45 : N_higher_{std::move(N1)}, N_lower_{std::move(N2)}, weight_{weight}
46 {
47 }
48
49 double weight() const { return weight_; }
50
51 Eigen::Ref<const Eigen::RowVectorXd> N(unsigned order) const
52 {
53 if (order < 2)
54 {
55 return N_lower_;
56 }
57
58 if (order == 2)
59 {
60 return N_higher_;
61 }
62
64 "Only shape functions up to order 2 are supported currently. You "
65 "have requested order {}. There might be an error in the OGS "
66 "project file.",
67 order);
68 }
69
70 ShapeMatrix const& NHigherOrder() const { return N_higher_; }
71
72private:
73 ShapeMatrix N_higher_;
74 LowerOrderShapeMatrix N_lower_;
75 double weight_;
76};
77
80template <typename ShapeFunction, typename ShapeMatrix>
81struct NsAndWeight<ShapeFunction, ShapeFunction, ShapeMatrix, ShapeMatrix>
82{
83 static_assert(ShapeFunction::ORDER < 2,
84 "For higher order shape functions the general case of this "
85 "class template should be used. See above.");
86
87 NsAndWeight(ShapeMatrix&& N, double const weight)
88 : N_{std::move(N)}, weight_{weight}
89 {
90 }
91
92 Eigen::Ref<const Eigen::RowVectorXd> N(unsigned order) const
93 {
94 if (order >= 2)
95 {
97 "Only shape functions of order < 2 are available in the "
98 "current setting. You have requested order {}. Maybe there is "
99 "an error in the OGS project file.",
100 order);
101 }
102
103 return N_;
104 }
105
106 ShapeMatrix const& NHigherOrder() const { return N_; }
107
108 double weight() const { return weight_; }
109
110private:
111 ShapeMatrix N_;
112 double weight_;
113};
114
116template <typename ShapeFunction, typename LowerOrderShapeFunction,
117 int GlobalDim>
119{
121 using ShapeMatrix = typename ShapeMatrixPolicy::ShapeMatrices::ShapeType;
122
126 typename LowerOrderShapeMatrixPolicy::ShapeMatrices::ShapeType;
127
130 ShapeFunction, LowerOrderShapeFunction, ShapeMatrix,
132};
133
141template <typename ShapeFunction, typename LowerOrderShapeFunction,
142 int GlobalDim, typename IntegrationMethod>
144 bool const is_axially_symmetric,
145 IntegrationMethod const& integration_method)
146{
147 using Traits =
149 using VecOfNsAndWeight = std::vector<typename Traits::NsAndWeight>;
150
151 VecOfNsAndWeight nss_and_weights;
152 nss_and_weights.reserve(integration_method.getNumberOfPoints());
153
154 auto sms =
155 NumLib::initShapeMatrices<ShapeFunction,
156 typename Traits::ShapeMatrixPolicy, GlobalDim,
158 element, is_axially_symmetric, integration_method);
159
160 if constexpr (std::is_same_v<ShapeFunction, LowerOrderShapeFunction>)
161 {
162 static_assert(ShapeFunction::ORDER < 2,
163 "We do not expect higher order shape functions here. "
164 "Something must have gone terribly wrong.");
165
166 for (unsigned ip = 0; ip < sms.size(); ++ip)
167 {
168 auto& sm = sms[ip];
169 double const w =
170 sm.detJ * sm.integralMeasure *
171 integration_method.getWeightedPoint(ip).getWeight();
172
173 nss_and_weights.emplace_back(std::move(sm.N), w);
174 }
175 }
176 else
177 {
178 auto sms_lower = NumLib::initShapeMatrices<
179 LowerOrderShapeFunction,
180 typename Traits::LowerOrderShapeMatrixPolicy, GlobalDim,
181 NumLib::ShapeMatrixType::N>(element, is_axially_symmetric,
182 integration_method);
183
184 for (unsigned ip = 0; ip < sms.size(); ++ip)
185 {
186 auto& sm = sms[ip];
187
188 // Note: we use det(J) of the higher order shape function. For
189 // linear geometries (i.e., higher order elements but with flat
190 // surfaces) there should be no difference.
191 double const w =
192 sm.detJ * sm.integralMeasure *
193 integration_method.getWeightedPoint(ip).getWeight();
194
195 nss_and_weights.emplace_back(std::move(sm.N),
196 std::move(sms_lower[ip].N), w);
197 }
198 }
199
200 return nss_and_weights;
201}
202
203} // namespace ProcessLib::Python
Definition of the Element class.
#define OGS_FATAL(...)
Definition Error.h:26
std::vector< typename ShapeMatricesType::ShapeMatrices, Eigen::aligned_allocator< typename ShapeMatricesType::ShapeMatrices > > initShapeMatrices(MeshLib::Element const &e, bool const is_axially_symmetric, IntegrationMethod const &integration_method)
@ N_J
calculates N, dNdr, J, and detJ
auto computeNsAndWeights(MeshLib::Element const &element, bool const is_axially_symmetric, IntegrationMethod const &integration_method)
NsAndWeight(ShapeMatrix &&N1, LowerOrderShapeMatrix &&N2, double const weight)
Definition NsAndWeight.h:43
Eigen::Ref< const Eigen::RowVectorXd > N(unsigned order) const
Definition NsAndWeight.h:51
Collects common type aliases needed when working with NsAndWeight.
typename LowerOrderShapeMatrixPolicy::ShapeMatrices::ShapeType LowerOrderShapeMatrix
typename ShapeMatrixPolicy::ShapeMatrices::ShapeType ShapeMatrix