FrictionQPotFEM 0.23.3
Loading...
Searching...
No Matches
TyingsPeriodic.h
Go to the documentation of this file.
1
9#ifndef GOOSEFEM_TYINGSPERIODIC_H
10#define GOOSEFEM_TYINGSPERIODIC_H
11
12#include "Mesh.h"
13#include "config.h"
14
15#include <Eigen/Eigen>
16#include <Eigen/Sparse>
17
18namespace GooseFEM {
19
23namespace Tyings {
24
46class Periodic {
47public:
48 Periodic() = default;
49
62 template <class C, class D, class S, class T>
63 Periodic(const C& coor, const D& dofs, const S& control_dofs, const T& nodal_tyings)
64 : Periodic(coor, dofs, control_dofs, nodal_tyings, xt::eval(xt::empty<size_t>({0})))
65 {
66 }
67
82 template <class C, class D, class S, class T, class U>
84 const C& coor,
85 const D& dofs,
86 const S& control_dofs,
87 const T& nodal_tyings,
88 const U& iip)
89 {
90 m_tyings = nodal_tyings;
91 m_coor = coor;
92 m_ndim = m_coor.shape(1);
93 m_nties = m_tyings.shape(0);
94
95 GOOSEFEM_ASSERT(xt::has_shape(m_tyings, {m_nties, size_t(2)}));
96 GOOSEFEM_ASSERT(xt::has_shape(control_dofs, {m_ndim, m_ndim}));
97 GOOSEFEM_ASSERT(xt::has_shape(dofs, m_coor.shape()));
98 GOOSEFEM_ASSERT(xt::amax(control_dofs)() <= xt::amax(dofs)());
99 GOOSEFEM_ASSERT(xt::amin(control_dofs)() >= xt::amin(dofs)());
100 GOOSEFEM_ASSERT(xt::amax(iip)() <= xt::amax(dofs)());
101 GOOSEFEM_ASSERT(xt::amin(iip)() >= xt::amin(dofs)());
102 GOOSEFEM_ASSERT(xt::amax(nodal_tyings)() < m_coor.shape(0));
103
104 array_type::tensor<size_t, 1> dependent = xt::view(m_tyings, xt::all(), 1);
105 array_type::tensor<size_t, 2> dependent_dofs =
106 xt::view(dofs, xt::keep(dependent), xt::all());
107 U iid = xt::flatten(dependent_dofs);
108 U iii = xt::setdiff1d(dofs, iid);
109 U iiu = xt::setdiff1d(iii, iip);
110
111 m_nnu = iiu.size();
112 m_nnp = iip.size();
113 m_nni = iii.size();
114 m_nnd = iid.size();
115
116 GooseFEM::Mesh::Reorder reorder({iiu, iip, iid});
117
118 m_dofs = reorder.apply(dofs);
119 m_control = reorder.apply(control_dofs);
120 }
121
125 size_t nnd() const
126 {
127 return m_nnd;
128 }
129
133 size_t nni() const
134 {
135 return m_nni;
136 }
137
141 size_t nnu() const
142 {
143 return m_nnu;
144 }
145
149 size_t nnp() const
150 {
151 return m_nnp;
152 }
153
158 {
159 return m_dofs;
160 }
161
167 {
168 return m_control;
169 }
170
179 {
180 return m_tyings;
181 }
182
189 {
190 return xt::arange<size_t>(m_nni, m_nni + m_nnd);
191 }
192
199 {
200 return xt::arange<size_t>(m_nni);
201 }
202
209 {
210 return xt::arange<size_t>(m_nnu);
211 }
212
219 {
220 return xt::arange<size_t>(m_nnp) + m_nnu;
221 }
222
237 Eigen::SparseMatrix<double> Cdi() const
238 {
239 std::vector<Eigen::Triplet<double>> data;
240
241 data.reserve(m_nties * m_ndim * (m_ndim + 1));
242
243 for (size_t i = 0; i < m_nties; ++i) {
244 for (size_t j = 0; j < m_ndim; ++j) {
245
246 size_t ni = m_tyings(i, 0);
247 size_t nd = m_tyings(i, 1);
248
249 data.push_back(Eigen::Triplet<double>(i * m_ndim + j, m_dofs(ni, j), +1.0));
250
251 for (size_t k = 0; k < m_ndim; ++k) {
252 data.push_back(Eigen::Triplet<double>(
253 i * m_ndim + j, m_control(j, k), m_coor(nd, k) - m_coor(ni, k)));
254 }
255 }
256 }
257
258 Eigen::SparseMatrix<double> Cdi;
259 Cdi.resize(m_nnd, m_nni);
260 Cdi.setFromTriplets(data.begin(), data.end());
261
262 return Cdi;
263 }
264
270 Eigen::SparseMatrix<double> Cdu() const
271 {
272 std::vector<Eigen::Triplet<double>> data;
273
274 data.reserve(m_nties * m_ndim * (m_ndim + 1));
275
276 for (size_t i = 0; i < m_nties; ++i) {
277 for (size_t j = 0; j < m_ndim; ++j) {
278
279 size_t ni = m_tyings(i, 0);
280 size_t nd = m_tyings(i, 1);
281
282 if (m_dofs(ni, j) < m_nnu) {
283 data.push_back(Eigen::Triplet<double>(i * m_ndim + j, m_dofs(ni, j), +1.0));
284 }
285
286 for (size_t k = 0; k < m_ndim; ++k) {
287 if (m_control(j, k) < m_nnu) {
288 data.push_back(Eigen::Triplet<double>(
289 i * m_ndim + j, m_control(j, k), m_coor(nd, k) - m_coor(ni, k)));
290 }
291 }
292 }
293 }
294
295 Eigen::SparseMatrix<double> Cdu;
296 Cdu.resize(m_nnd, m_nnu);
297 Cdu.setFromTriplets(data.begin(), data.end());
298
299 return Cdu;
300 }
301
307 Eigen::SparseMatrix<double> Cdp() const
308 {
309 std::vector<Eigen::Triplet<double>> data;
310
311 data.reserve(m_nties * m_ndim * (m_ndim + 1));
312
313 for (size_t i = 0; i < m_nties; ++i) {
314 for (size_t j = 0; j < m_ndim; ++j) {
315
316 size_t ni = m_tyings(i, 0);
317 size_t nd = m_tyings(i, 1);
318
319 if (m_dofs(ni, j) >= m_nnu) {
320 data.push_back(
321 Eigen::Triplet<double>(i * m_ndim + j, m_dofs(ni, j) - m_nnu, +1.0));
322 }
323
324 for (size_t k = 0; k < m_ndim; ++k) {
325 if (m_control(j, k) >= m_nnu) {
326 data.push_back(Eigen::Triplet<double>(
327 i * m_ndim + j,
328 m_control(j, k) - m_nnu,
329 m_coor(nd, k) - m_coor(ni, k)));
330 }
331 }
332 }
333 }
334
335 Eigen::SparseMatrix<double> Cdp;
336 Cdp.resize(m_nnd, m_nnp);
337 Cdp.setFromTriplets(data.begin(), data.end());
338
339 return Cdp;
340 }
341
342private:
343 size_t m_nnu;
344 size_t m_nnp;
345 size_t m_nni;
346 size_t m_nnd;
347 size_t m_ndim;
348 size_t m_nties;
353};
354
358class Control {
359public:
360 Control() = default;
361
370 template <class C, class N>
371 Control(const C& coor, const N& dofs)
372 {
373 GOOSEFEM_ASSERT(coor.shape().size() == 2);
374 GOOSEFEM_ASSERT(coor.shape() == dofs.shape());
375
376 m_coor = coor;
377 m_dofs = dofs;
378
379 size_t nnode = coor.shape(0);
380 size_t ndim = coor.shape(1);
381
382 m_control_dofs = xt::arange<size_t>(ndim * ndim).reshape({ndim, ndim});
383 m_control_dofs += xt::amax(dofs)() + 1;
384
385 m_control_nodes = nnode + xt::arange<size_t>(ndim);
386
387 m_coor = xt::concatenate(xt::xtuple(coor, xt::zeros<double>({ndim, ndim})));
388 m_dofs = xt::concatenate(xt::xtuple(dofs, m_control_dofs));
389 }
390
397 {
398 return m_coor;
399 }
400
407 {
408 return m_dofs;
409 }
410
417 {
418 return m_control_dofs;
419 }
420
427 {
428 return m_control_nodes;
429 }
430
431private:
434 array_type::tensor<size_t, 2> m_control_dofs;
435 array_type::tensor<size_t, 1> m_control_nodes;
436};
437
438} // namespace Tyings
439} // namespace GooseFEM
440
441#endif
Generic mesh operations.
Reorder to lowest possible index, in specific order.
Definition: Mesh.h:2387
Add control nodes to an existing system.
const array_type::tensor< double, 2 > & coor() const
Nodal coordinates, for the system with control nodes added to it.
const array_type::tensor< size_t, 2 > & dofs() const
DOF-numbers per node, for the system with control nodes added to it.
const array_type::tensor< size_t, 1 > & controlNodes() const
Node-numbers of the control nodes.
Control(const C &coor, const N &dofs)
Constructor.
const array_type::tensor< size_t, 2 > & controlDofs() const
DOF-numbers of each control node.
Nodal tyings per periodic boundary conditions.
Periodic(const C &coor, const D &dofs, const S &control_dofs, const T &nodal_tyings)
Constructor.
Periodic(const C &coor, const D &dofs, const S &control_dofs, const T &nodal_tyings, const U &iip)
Constructor.
Eigen::SparseMatrix< double > Cdp() const
Prescribed part of the partitioned tying matrix, see Cdi().
Eigen::SparseMatrix< double > Cdu() const
Unknown part of the partitioned tying matrix, see Cdi().
array_type::tensor< size_t, 1 > iip() const
Independent prescribed DOFs.
array_type::tensor< size_t, 1 > iiu() const
Independent unknown DOFs.
array_type::tensor< size_t, 1 > iid() const
Dependent DOFs.
array_type::tensor< size_t, 1 > iii() const
Independent DOFs.
const array_type::tensor< size_t, 2 > & dofs() const
const array_type::tensor< size_t, 2 > & control() const
Eigen::SparseMatrix< double > Cdi() const
Return tying matrix such as to get the dependent DOFs from the independent DOFs as follows.
const array_type::tensor< size_t, 2 > & nodal_tyings() const
Return the applied nodal tyings.
#define GOOSEFEM_ASSERT(expr)
All assertions are implementation as::
Definition: config.h:96
xt::xtensor< T, N > tensor
Fixed (static) rank array.
Definition: config.h:176
Toolbox to perform finite element computations.
Definition: Allocate.h:14