HArD::Core2D
Hybrid Arbitrary Degree::Core 2D - Library to implement 2D schemes with edge and cell polynomials as unknowns
Loading...
Searching...
No Matches
xdivdiv-test.hpp
Go to the documentation of this file.
1// Test of the xdivdiv class
2//
3// Author: Daniele Di Pietro (daniele.di-pietro@umontpellier.fr), Jerome Droniou (jerome.droniou@monash.edu)
4//
5
6
7#ifndef XDIVDIV_TEST_HPP
8#define XDIVDIV_TEST_HPP
9
10#include <iostream>
11
12#include <boost/math/constants/constants.hpp>
13
14#include <Eigen/Sparse>
15#include <unsupported/Eigen/SparseExtra>
16
17#include <mesh.hpp>
18
19#include <xdivdiv.hpp>
20
21
22namespace HArDCore2D
23{
24
25
27 struct KirchhoffLove
28 {
29 typedef Eigen::SparseMatrix<double> SystemMatrixType;
30
31 typedef std::function<double(const Eigen::Vector2d &)> ForcingTermType;
32 typedef std::function<double(const Eigen::Vector2d &)> DeflectionType;
33 typedef std::function<Eigen::Matrix2d(const Eigen::Vector2d &)> MomentTensorType;
34 typedef std::function<double(const Eigen::Vector2d &, const Edge &)> MomentTensorEdgeDerivativeType;
35
36
38 KirchhoffLove(
39 const PlatesCore & platescore,
40 bool use_threads,
41 std::ostream & output = std::cout
42 );
43
45 inline size_t dimensionSpace() const
46 {
47 return m_xdivdiv.dimension() + m_platescore.mesh().n_cells() * PolynomialSpaceDimension<Cell>::Poly(m_platescore.degree()-2);
48 }
49
51 inline const XDivDiv & xDivDiv() const {
52 return m_xdivdiv;
53 }
54
56 inline const SystemMatrixType & systemMatrix() const {
57 return m_A;
58 }
59
61 inline SystemMatrixType & systemMatrix() {
62 return m_A;
63 }
64
66 inline const Eigen::VectorXd & systemVector() const {
67 return m_b;
68 }
69
71 inline Eigen::VectorXd & systemVector() {
72 return m_b;
73 }
74
76 inline const double & stabilizationParameter() const {
77 return m_stab_par;
78 }
79
81 inline double & stabilizationParameter() {
82 return m_stab_par;
83 }
84
85 private:
86 const PlatesCore & m_platescore;
87// bool m_use_threads;
88 std::ostream & m_output;
89
90 XDivDiv m_xdivdiv;
91
92 SystemMatrixType m_A;
93 Eigen::VectorXd m_b;
94
95 double m_stab_par;
96 };
97
98 //------------------------------------------------------------------------------
99 // Exact solutions
100 //------------------------------------------------------------------------------
101
102 static const double PI = boost::math::constants::pi<double>();
103 using std::sin;
104 using std::cos;
105 using std::exp;
106 using std::pow;
107
108 //------------------------------------------------------------------------------
109 // Trigonometric
110
111 static KirchhoffLove::DeflectionType
112 trigonometric_u = [](const VectorRd & x) -> double {
113 return sin(PI*x(0))*sin(PI*x(1));
114 };
115
116 static KirchhoffLove::ForcingTermType
117 trigonometric_f = [](const VectorRd & x) -> double {
118 return 4.*pow(PI, 4)*sin(PI*x(0))*sin(PI*x(1));
119 };
120
121 static KirchhoffLove::MomentTensorType
122 trigonometric_sigma = [](const VectorRd & x) -> Eigen::Matrix2d {
123 Eigen::Matrix2d M;
124 M.row(0) << -sin(PI*x(0))*sin(PI*x(1)), cos(PI*x(0))*cos(PI*x(1));
125 M.row(1) << cos(PI*x(0))*cos(PI*x(1)), -sin(PI*x(0))*sin(PI*x(1));
126
127 return -pow(PI, 2) * M;
128 };
129
130 static KirchhoffLove::MomentTensorEdgeDerivativeType
131 trigonometric_sigma_DE = [](const Eigen::Vector2d & x, const Edge & E) -> double {
132
133 Eigen::Vector2d grad_sigma11 {
134 -cos(PI*x(0))*sin(PI*x(1)), -sin(PI*x(0))*cos(PI*x(1))
135 };
136 grad_sigma11 *= -pow(PI, 3);
137
138 Eigen::Vector2d grad_sigma12 {
139 -sin(PI*x(0))*cos(PI*x(1)), -sin(PI*x(1))*cos(PI*x(0))
140 };
141 grad_sigma12 *= -pow(PI, 3);
142
143 Eigen::Vector2d grad_sigma22 {
144 -cos(PI*x(0))*sin(PI*x(1)), -sin(PI*x(0))*cos(PI*x(1))
145 };
146 grad_sigma22 *= -pow(PI, 3);
147
148 Eigen::Vector2d tE = E.tangent();
149 Eigen::Vector2d nE = E.normal();
150
151 Eigen::Matrix2d dtE_sigma;
152 dtE_sigma.row(0) << grad_sigma11.dot(tE), grad_sigma12.dot(tE);
153 dtE_sigma.row(1) << grad_sigma12.dot(tE), grad_sigma22.dot(tE);
154
155 Eigen::Vector2d div_sigma {
156 grad_sigma11(0) + grad_sigma12(1),
157 grad_sigma12(0) + grad_sigma22(1)
158 };
159
160 return (dtE_sigma*nE).dot(tE) + div_sigma.dot(nE);
161 };
162
163 //------------------------------------------------------------------------------
164 // Simple
165
166 static KirchhoffLove::DeflectionType
167 simple_u = [](const VectorRd & x) -> double {
168 return 1.;
169 };
170
171// static KirchhoffLove::GradientDeflectionType
172// simple_grad_u = [](const VectorRd & x) -> VectorRd {
173// return VectorRd(0., 0.);
174// };
175
176 static KirchhoffLove::ForcingTermType
177 simple_f = [](const VectorRd & x) -> double {
178 return 0.;
179 };
180
181 static KirchhoffLove::MomentTensorType
182 simple_sigma = [](const VectorRd & x) -> Eigen::Matrix2d {
183 // Here, M=hess(u) (modulo factor pi^2) and sigma is the opposite
184 Eigen::Matrix2d M;
185 M.row(0) << 0., 0.;
186 M.row(1) << 0., 0.;
187
188 return M;
189 };
190
191 static KirchhoffLove::MomentTensorEdgeDerivativeType
192 simple_sigma_DE = [](const Eigen::Vector2d & x, const Edge & E) -> double {
193
194 Eigen::Vector2d grad_sigma11 {
195 0., 0.
196 };
197 grad_sigma11 *= -pow(PI, 3);
198
199 Eigen::Vector2d grad_sigma12 {
200 0., 0.
201 };
202 grad_sigma12 *= -pow(PI, 3);
203
204 Eigen::Vector2d grad_sigma22 {
205 0., 0.
206 };
207 grad_sigma22 *= -pow(PI, 3);
208
209 Eigen::Vector2d tE = E.tangent();
210 Eigen::Vector2d nE = E.normal();
211
212 Eigen::Matrix2d dtE_sigma;
213 dtE_sigma.row(0) << grad_sigma11.dot(tE), grad_sigma12.dot(tE);
214 dtE_sigma.row(1) << grad_sigma12.dot(tE), grad_sigma22.dot(tE);
215
216 Eigen::Vector2d div_sigma {
217 grad_sigma11(0) + grad_sigma12(1),
218 grad_sigma12(0) + grad_sigma22(1)
219 };
220
221 return (dtE_sigma*nE).dot(tE) + div_sigma.dot(nE);
222 };
223
224 //------------------------------------------------------------------------------
225 // Quartic
226
227 static KirchhoffLove::DeflectionType
228 quartic_u = [](const VectorRd & x) -> double {
229 return pow(x(0),2)*pow(1.-x(0),2) + pow(x(1),2)*pow(1.-x(1),2);
230 };
231
232 static KirchhoffLove::ForcingTermType
233 quartic_f = [](const VectorRd & x) -> double {
234 return 48.;
235 };
236
237 static KirchhoffLove::MomentTensorType
238 quartic_sigma = [](const VectorRd & x) -> Eigen::Matrix2d {
239
240 Eigen::Matrix2d M;
241 M <<
242 2.*(pow(x(0),2)+4.*x(0)*(x(0)-1.)+pow(1.-x(0),2)), 0.,
243 0., 2.*(pow(x(1),2)+4.*x(1)*(x(1)-1.)+pow(1.-x(1),2));
244
245 return -M;
246 };
247
248 static KirchhoffLove::MomentTensorEdgeDerivativeType
249 quartic_sigma_DE = [](const Eigen::Vector2d & x, const Edge & E) -> double {
250
251 Eigen::Vector2d grad_sigma11 { -24.*x(0)+12., 0. };
252 Eigen::Vector2d grad_sigma12 { 0., 0. };
253 Eigen::Vector2d grad_sigma22 { 0., -24.*x(1)+12. };
254
255 Eigen::Vector2d tE = E.tangent();
256 Eigen::Vector2d nE = E.normal();
257
258 Eigen::Matrix2d dtE_sigma;
259 dtE_sigma <<
260 grad_sigma11.dot(tE), grad_sigma12.dot(tE),
261 grad_sigma12.dot(tE), grad_sigma22.dot(tE);
262
263 Eigen::Vector2d div_sigma {
264 grad_sigma11(0) + grad_sigma12(1),
265 grad_sigma12(0) + grad_sigma22(1)
266 };
267
268 return (dtE_sigma*nE).dot(tE) + div_sigma.dot(nE);
269 };
270
271 //------------------------------------------------------------------------------
272 // Biquartic
273
274 static KirchhoffLove::DeflectionType
275 biquartic_u = [](const VectorRd & x) -> double {
276 return pow(x(0),2)*pow(1.-x(0),2)*pow(x(1),2)*pow(1.-x(1),2);
277 };
278
279 static KirchhoffLove::ForcingTermType
280 biquartic_f = [](const VectorRd & x) -> double {
281 return 24*pow(x(0), 2)*pow(x(0) - 1, 2) + 32*x(0)*x(1)*(x(0) - 1)*(x(1) - 1) + 8*x(0)*x(1)*(x(0) - 1)*(4*x(1) - 2) + 8*x(0)*x(1)*(4*x(0) - 2)*(x(1) - 1) + 8*x(0)*x(1)*(x(0)*x(1) + x(0)*(x(1) - 1) + x(1)*(x(0) - 1) + (x(0) - 1)*(x(1) - 1)) + 8*x(0)*(x(0) - 1)*(x(1) - 1)*(4*x(1) - 2) + 8*x(0)*(x(1) - 1)*(x(0)*x(1) + x(0)*(x(1) - 1) + x(1)*(x(0) - 1) + (x(0) - 1)*(x(1) - 1)) + 24*pow(x(1), 2)*pow(x(1) - 1, 2) + 8*x(1)*(x(0) - 1)*(4*x(0) - 2)*(x(1) - 1) + 8*x(1)*(x(0) - 1)*(x(0)*x(1) + x(0)*(x(1) - 1) + x(1)*(x(0) - 1) + (x(0) - 1)*(x(1) - 1)) + 2*(4*x(0) - 4)*(x(1) - 1)*(x(0)*x(1) + x(0)*(x(1) - 1) + x(1)*(x(0) - 1) + (x(0) - 1)*(x(1) - 1));
282 };
283
284 static KirchhoffLove::MomentTensorType
285 biquartic_sigma = [](const VectorRd & x) -> Eigen::Matrix2d {
286 Eigen::Matrix2d M;
287 M.row(0) << 2*pow(x(1), 2)*pow(x(1) - 1, 2)*(pow(x(0), 2) + 4*x(0)*(x(0) - 1) + pow(x(0) - 1, 2)), 4*x(0)*x(1)*(x(0) - 1)*(x(1) - 1)*(x(0)*x(1) + x(0)*(x(1) - 1) + x(1)*(x(0) - 1) + (x(0) - 1)*(x(1) - 1));
288 M.row(1) << 4*x(0)*x(1)*(x(0) - 1)*(x(1) - 1)*(x(0)*x(1) + x(0)*(x(1) - 1) + x(1)*(x(0) - 1) + (x(0) - 1)*(x(1) - 1)), 2*pow(x(0), 2)*pow(x(0) - 1, 2)*(pow(x(1), 2) + 4*x(1)*(x(1) - 1) + pow(x(1) - 1, 2));
289 return -M;
290 };
291
292 static KirchhoffLove::MomentTensorEdgeDerivativeType
293 biquartic_sigma_DE = [](const Eigen::Vector2d & x, const Edge & E) -> double {
294
295 Eigen::Vector2d grad_sigma11 {
296 2*pow(x(1), 2)*(12*x(0) - 6)*pow(x(1) - 1, 2),
297 2*pow(x(1), 2)*(2*x(1) - 2)*(pow(x(0), 2) + 4*x(0)*(x(0) - 1) + pow(x(0) - 1, 2)) + 4*x(1)*pow(x(1) - 1, 2)*(pow(x(0), 2) + 4*x(0)*(x(0) - 1) + pow(x(0) - 1, 2))
298 };
299
300 Eigen::Vector2d grad_sigma12 {
301 4*x(0)*x(1)*(x(0) - 1)*(x(1) - 1)*(4*x(1) - 2) + 4*x(0)*x(1)*(x(1) - 1)*(x(0)*x(1) + x(0)*(x(1) - 1) + x(1)*(x(0) - 1) + (x(0) - 1)*(x(1) - 1)) + 4*x(1)*(x(0) - 1)*(x(1) - 1)*(x(0)*x(1) + x(0)*(x(1) - 1) + x(1)*(x(0) - 1) + (x(0) - 1)*(x(1) - 1)),
302 4*x(0)*x(1)*(x(0) - 1)*(4*x(0) - 2)*(x(1) - 1) + 4*x(0)*x(1)*(x(0) - 1)*(x(0)*x(1) + x(0)*(x(1) - 1) + x(1)*(x(0) - 1) + (x(0) - 1)*(x(1) - 1)) + 4*x(0)*(x(0) - 1)*(x(1) - 1)*(x(0)*x(1) + x(0)*(x(1) - 1) + x(1)*(x(0) - 1) + (x(0) - 1)*(x(1) - 1))
303 };
304
305 Eigen::Vector2d grad_sigma22 {
306 2*pow(x(0), 2)*(2*x(0) - 2)*(pow(x(1), 2) + 4*x(1)*(x(1) - 1) + pow(x(1) - 1, 2)) + 4*x(0)*pow(x(0) - 1, 2)*(pow(x(1), 2) + 4*x(1)*(x(1) - 1) + pow(x(1) - 1, 2)),
307 2*pow(x(0), 2)*pow(x(0) - 1, 2)*(12*x(1) - 6)
308 };
309
310 // The previous calculations are actually grad_hess(u)
311 grad_sigma11 *= -1.;
312 grad_sigma12 *= -1.;
313 grad_sigma22 *= -1.;
314
315 Eigen::Vector2d tE = E.tangent();
316 Eigen::Vector2d nE = E.normal();
317
318 Eigen::Matrix2d dtE_sigma;
319 dtE_sigma.row(0) << grad_sigma11.dot(tE), grad_sigma12.dot(tE);
320 dtE_sigma.row(1) << grad_sigma12.dot(tE), grad_sigma22.dot(tE);
321
322 Eigen::Vector2d div_sigma {
323 grad_sigma11(0) + grad_sigma12(1),
324 grad_sigma12(0) + grad_sigma22(1)
325 };
326
327 return (dtE_sigma*nE).dot(tE) + div_sigma.dot(nE);
328 };
329
330} // end of namespace HArDCore2D
331
332#endif
HArDCore2D::PlatesCore
Construct all polynomial spaces for the plates sequence.
Definition platescore.hpp:25
HArDCore2D::XDivDiv
Discrete Hdivdiv space.
Definition xdivdiv.hpp:20
x
Create grid points x
Definition generate_cartesian_mesh.m:22
HArDCore2D::KirchhoffLove::SystemMatrixType
Eigen::SparseMatrix< double > SystemMatrixType
Definition bgg-klplate.hpp:54
HArDCore2D::PI
static const double PI
Definition bgg-klplate.hpp:214
HArDCore2D::KirchhoffLove::ForcingTermType
std::function< double(const Eigen::Vector2d &)> ForcingTermType
Definition bgg-klplate.hpp:55
HArDCore2D::VectorRd
Eigen::Vector2d VectorRd
Definition basis.hpp:55
HArDCore2D::LocalDOFSpace::dimension
size_t dimension() const
Returns the dimension of the global space (all DOFs for all geometric entities)
Definition localdofspace.hpp:61
HArDCore2D::trigonometric_u
static KirchhoffLove::DeflectionType trigonometric_u
Definition ddr-klplate.hpp:197
HArDCore2D::biquartic_u
static KirchhoffLove::DeflectionType biquartic_u
Definition ddr-klplate.hpp:346
HArDCore2D::KirchhoffLove::MomentTensorType
std::function< Eigen::Matrix2d(const Eigen::Vector2d &)> MomentTensorType
Definition ddr-klplate.hpp:42
HArDCore2D::quartic_u
static KirchhoffLove::DeflectionType quartic_u
Definition ddr-klplate.hpp:270
HArDCore2D::KirchhoffLove::MomentTensorEdgeDerivativeType
std::function< double(const Eigen::Vector2d &, const Edge &)> MomentTensorEdgeDerivativeType
Definition ddr-klplate.hpp:43
HArDCore2D::KirchhoffLove::DeflectionType
std::function< double(const Eigen::Vector2d &)> DeflectionType
Definition ddr-klplate.hpp:40
use_threads
bool use_threads
Definition HHO_DiffAdvecReac.hpp:45
Mesh2D::Mesh::n_cells
std::size_t n_cells() const
number of cells in the mesh.
Definition Mesh2D.hpp:63
HArDCore2D::PlatesCore::mesh
const Mesh & mesh() const
Return a const reference to the mesh.
Definition platescore.hpp:66
HArDCore2D::PlatesCore::degree
const size_t & degree() const
Return the polynomial degree.
Definition platescore.hpp:72
array
depending on the Matrix Market format indicated by or array(dense array storage). The data will be duplicated % as appropriate if symmetry is indicated in the header. % % Optionally
HArDCore2D
Definition mhd-solutions.hpp:9
HArDCore2D::trigonometric_sigma
static KirchhoffLove::MomentTensorType trigonometric_sigma
Definition xdivdiv-test.hpp:122
HArDCore2D::biquartic_sigma_DE
static KirchhoffLove::MomentTensorEdgeDerivativeType biquartic_sigma_DE
Definition xdivdiv-test.hpp:293
HArDCore2D::quartic_sigma_DE
static KirchhoffLove::MomentTensorEdgeDerivativeType quartic_sigma_DE
Definition xdivdiv-test.hpp:249
HArDCore2D::simple_u
static KirchhoffLove::DeflectionType simple_u
Definition xdivdiv-test.hpp:167
HArDCore2D::trigonometric_f
static QuadRot::ForcingTermType trigonometric_f
Definition ddr-quadrot.hpp:367
HArDCore2D::simple_f
static KirchhoffLove::ForcingTermType simple_f
Definition xdivdiv-test.hpp:177
HArDCore2D::quartic_f
static QuadRot::ForcingTermType quartic_f
Definition ddr-quadrot.hpp:327
HArDCore2D::simple_sigma
static KirchhoffLove::MomentTensorType simple_sigma
Definition xdivdiv-test.hpp:182
HArDCore2D::trigonometric_sigma_DE
static KirchhoffLove::MomentTensorEdgeDerivativeType trigonometric_sigma_DE
Definition xdivdiv-test.hpp:131
HArDCore2D::simple_sigma_DE
static KirchhoffLove::MomentTensorEdgeDerivativeType simple_sigma_DE
Definition xdivdiv-test.hpp:192
HArDCore2D::quartic_sigma
static KirchhoffLove::MomentTensorType quartic_sigma
Definition xdivdiv-test.hpp:238
HArDCore2D::biquartic_f
static KirchhoffLove::ForcingTermType biquartic_f
Definition xdivdiv-test.hpp:280
HArDCore2D::biquartic_sigma
static KirchhoffLove::MomentTensorType biquartic_sigma
Definition xdivdiv-test.hpp:285
HArDCore2D::KirchhoffLove
Assemble a RM problem.
Definition bgg-klplate.hpp:53
HArDCore2D::KirchhoffLove::systemVector
Eigen::VectorXd & systemVector()
Returns the linear system right-hand side vector.
Definition xdivdiv-test.hpp:71
HArDCore2D::KirchhoffLove::stabilizationParameter
double & stabilizationParameter()
Returns the stabilization parameter.
Definition xdivdiv-test.hpp:81
HArDCore2D::KirchhoffLove::SystemMatrixType
Eigen::SparseMatrix< double > SystemMatrixType
Definition xdivdiv-test.hpp:29
HArDCore2D::KirchhoffLove::ForcingTermType
std::function< double(const Eigen::Vector2d &)> ForcingTermType
Definition xdivdiv-test.hpp:31
HArDCore2D::KirchhoffLove::systemMatrix
const SystemMatrixType & systemMatrix() const
Returns the linear system matrix.
Definition xdivdiv-test.hpp:56
HArDCore2D::KirchhoffLove::systemVector
const Eigen::VectorXd & systemVector() const
Returns the linear system right-hand side vector.
Definition xdivdiv-test.hpp:66
HArDCore2D::KirchhoffLove::systemMatrix
SystemMatrixType & systemMatrix()
Returns the linear system matrix.
Definition xdivdiv-test.hpp:61
HArDCore2D::KirchhoffLove::MomentTensorType
std::function< Eigen::Matrix2d(const Eigen::Vector2d &)> MomentTensorType
Definition xdivdiv-test.hpp:33
HArDCore2D::KirchhoffLove::xDivDiv
const XDivDiv & xDivDiv() const
Returns the space XDivDiv.
Definition xdivdiv-test.hpp:51
HArDCore2D::KirchhoffLove::MomentTensorEdgeDerivativeType
std::function< double(const Eigen::Vector2d &, const Edge &)> MomentTensorEdgeDerivativeType
Definition xdivdiv-test.hpp:34
HArDCore2D::KirchhoffLove::stabilizationParameter
const double & stabilizationParameter() const
Returns the stabilization parameter.
Definition xdivdiv-test.hpp:76
HArDCore2D::KirchhoffLove::DeflectionType
std::function< double(const Eigen::Vector2d &)> DeflectionType
Definition xdivdiv-test.hpp:32
HArDCore2D::KirchhoffLove::dimensionSpace
size_t dimensionSpace() const
Returns the dimension of the rotation + displacement space (with BC)
Definition xdivdiv-test.hpp:45
HArDCore2D::PolynomialSpaceDimension
Basis dimensions for various polynomial spaces on edges/faces/elements (when relevant): Pk,...
Definition polynomialspacedimension.hpp:55
Generated by doxygen 1.9.8