G+Smo  25.01.0
Geometry + Simulation Modules
 
Loading...
Searching...
No Matches
gsALE.hpp
Go to the documentation of this file.
1
15#pragma once
16
17#include <gsElasticity/gsALE.h>
18
19#include <gsElasticity/gsElasticityAssembler.h>
20#include <gsElasticity/gsElPoissonAssembler.h>
21#include <gsElasticity/gsBiharmonicAssembler.h>
22#include <gsElasticity/gsIterative.h>
23#include <gsElasticity/gsGeoUtils.h>
24#include <gsCore/gsConstantFunction.h>
25
26namespace gismo
27{
28
29template <class T>
30gsALE<T>::gsALE(gsMultiPatch<T> & geometry, const gsMultiPatch<T> & displacement,
31 const gsBoundaryInterface & interfaceS2M, ale_method::method method)
32 : disp(displacement),
33 m_interface(interfaceS2M),
34 methodALE(method),
35 m_options(defaultOptions()),
36 initialized(false),
37 hasSavedState(false)
38{
39 // create input for the assembler
40 gsMultiBasis<T> basis(geometry);
41 gsBoundaryConditions<T> bcInfo;
42 for (gsMultiPatch<>::const_biterator it = geometry.bBegin(); it != geometry.bEnd(); ++it)
43 for (index_t d = 0; d < geometry.parDim(); ++d)
44 bcInfo.addCondition(it->patch,it->side(),condition_type::dirichlet,0,d);
45 gsConstantFunction<T> rhs(gsVector<>::Zero(geometry.parDim()),geometry.parDim());
46
47 // define assembler according to the method
48 if (methodALE == ale_method::TINE || methodALE == ale_method::TINE_StVK || methodALE == ale_method::ILE || methodALE == ale_method::LE)
49 {
50 assembler = typename gsBaseAssembler<T>::uPtr(new gsElasticityAssembler<T>(geometry,basis,bcInfo,rhs));
51 if (methodALE == ale_method::TINE || methodALE == ale_method::TINE_StVK)
52 {
53 assembler->options().setSwitch("Check",false);
54 solverNL = typename gsIterative<T>::uPtr(new gsIterative<T>(*assembler,
55 gsMatrix<>::Zero(assembler->numDofs(),1),
56 assembler->allFixedDofs()));
57 solverNL->options().setInt("Verbosity",solver_verbosity::none);
58 solverNL->options().setInt("Solver",linear_solver::LDLT);
59 }
60 if (methodALE == ale_method::TINE)
61 assembler->options().setInt("MaterialLaw",material_law::neo_hooke_ln);
62 if (methodALE == ale_method::TINE_StVK)
63 assembler->options().setInt("MaterialLaw",material_law::saint_venant_kirchhoff);
64
65 assembler->constructSolution(gsMatrix<T>::Zero(assembler->numDofs(),1),assembler->allFixedDofs(),ALEdisp);
66 }
67 else if (methodALE == ale_method::HE || methodALE == ale_method::IHE)
68 {
69 assembler = typename gsBaseAssembler<T>::uPtr(new gsElPoissonAssembler<T>(geometry,basis,bcInfo,rhs));
70 assembler->constructSolution(gsMatrix<T>::Zero(assembler->numDofs(),geometry.parDim()),assembler->allFixedDofs(),ALEdisp);
71 }
72 else if (methodALE == ale_method::BHE || methodALE == ale_method::IBHE)
73 {
74 gsBoundaryConditions<T> bcInfoBHE;
75 for (gsMultiPatch<>::const_biterator it = geometry.bBegin(); it != geometry.bEnd(); ++it)
76 bcInfoBHE.addCondition(it->patch,it->side(),condition_type::dirichlet,0);
77 assembler = typename gsBaseAssembler<T>::uPtr(new gsBiharmonicAssembler<T>(geometry,basis,bcInfoBHE,rhs));
78 assembler->constructSolution(gsMatrix<T>::Zero(assembler->numDofs(),geometry.parDim()),assembler->allFixedDofs(),ALEdisp);
79 }
80 else
81 for (size_t p = 0; p < geometry.nPatches(); ++p)
82 {
83 ALEdisp.addPatch(geometry.patch(p).clone());
84 ALEdisp.patch(p).coefs() *= 0.;
85 }
86}
87
88template <class T>
89gsOptionList gsALE<T>::defaultOptions()
90{
91 gsOptionList opt;
92 opt.addReal("PoissonsRatio","Poisson's ratio of the material (only for elasticity-based methods)",0.4);
93 opt.addReal("LocalStiff","Stiffening degree for the Jacobian-based local stiffening",0.);
94 opt.addSwitch("Check","Check bijectivity of the resulting ALE displacement field",true);
95 opt.addInt("NumIter","Number of iterations for nonlinear methods",1);
96 return opt;
97}
98
99
100template <class T>
101void gsALE<T>::constructSolution(gsMultiPatch<T> & solution) const
102{
103 solution.clear();
104 for (size_t p = 0; p < ALEdisp.nPatches(); ++p)
105 solution.addPatch(ALEdisp.patch(p).clone());
106}
107
108template <class T>
109void gsALE<T>::initialize()
110{
111 assembler->options().setReal("LocalStiff",m_options.getReal("LocalStiff"));
112 if (methodALE == ale_method::LE || methodALE == ale_method::ILE || methodALE == ale_method::TINE || methodALE == ale_method::TINE_StVK)
113 assembler->options().setReal("PoissonsRatio",m_options.getReal("PoissonsRatio"));
114 if (methodALE == ale_method::LE || methodALE == ale_method::HE || methodALE == ale_method::BHE)
115 assembler->assemble(true);
116 if (methodALE == ale_method::TINE || methodALE == ale_method::TINE_StVK)
117 solverNL->options().setInt("MaxIters",m_options.getInt("NumIter"));
118
119 initialized = true;
120}
121
122template <class T>
123index_t gsALE<T>::updateMesh()
124{
125 if (!initialized)
126 initialize();
127
128 switch (methodALE)
129 {
130 case ale_method::HE: return linearMethod(); break;
131 case ale_method::IHE: return linearIncrementalMethod(); break;
132 case ale_method::LE: return linearMethod(); break;
133 case ale_method::ILE: return linearIncrementalMethod(); break;
134 case ale_method::TINE: return nonlinearMethod(); break;
135 case ale_method::TINE_StVK: return nonlinearMethod(); break;
136 case ale_method::BHE: return linearMethod(); break;
137 case ale_method::IBHE: return linearIncrementalMethod(); break;
138 default: return -1;
139 }
140}
141
142template <class T>
143index_t gsALE<T>::linearMethod()
144{
145
146 for (size_t i = 0; i < m_interface.sidesA.size(); ++i)
147 assembler->setFixedDofs(m_interface.sidesB[i].patch,
148 m_interface.sidesB[i].side(),
149 disp.patch(m_interface.sidesA[i].patch).boundary(m_interface.sidesA[i].side())->coefs(),
150 methodALE == ale_method::LE ? false : true);
151 assembler->eliminateFixedDofs();
152
153#ifdef GISMO_WITH_PARDISO
154 gsSparseSolver<>::PardisoLDLT solver(assembler->matrix());
155 gsMatrix<> solVector = solver.solve(assembler->rhs());
156#else
157 gsSparseSolver<>::SimplicialLDLT solver(assembler->matrix());
158 gsMatrix<> solVector = solver.solve(assembler->rhs());
159#endif
160
161 assembler->constructSolution(solVector,assembler->allFixedDofs(),ALEdisp);
162 if (m_options.getSwitch("Check"))
163 return checkDisplacement(assembler->patches(),ALEdisp);
164 else
165 return -1;
166}
167
168template <class T>
169index_t gsALE<T>::linearIncrementalMethod()
170{
171 for (size_t i = 0; i < m_interface.sidesA.size(); ++i)
172 assembler->setFixedDofs(m_interface.sidesB[i].patch,
173 m_interface.sidesB[i].side(),
174 disp.patch(m_interface.sidesA[i].patch).boundary(m_interface.sidesA[i].side())->coefs() -
175 ALEdisp.patch(m_interface.sidesB[i].patch).boundary(m_interface.sidesB[i].side())->coefs(),
176 methodALE == ale_method::ILE ? false : true);
177 assembler->assemble();
178
179#ifdef GISMO_WITH_PARDISO
180 gsSparseSolver<>::PardisoLDLT solver(assembler->matrix());
181 gsMatrix<> solVector = solver.solve(assembler->rhs());
182#else
183 gsSparseSolver<>::SimplicialLDLT solver(assembler->matrix());
184 gsMatrix<> solVector = solver.solve(assembler->rhs());
185#endif
186
187 gsMultiPatch<T> ALEupdate;
188 assembler->constructSolution(solVector,assembler->allFixedDofs(),ALEupdate);
189 for (size_t p = 0; p < ALEupdate.nPatches(); ++p)
190 {
191 ALEdisp.patch(p).coefs() += ALEupdate.patch(p).coefs();
192 assembler->patches().patch(p).coefs() += ALEupdate.patch(p).coefs();
193 }
194 if (m_options.getSwitch("Check"))
195 return checkGeometry(assembler->patches());
196 else
197 return -1;
198}
199
200template <class T>
201index_t gsALE<T>::nonlinearMethod()
202{
203 for (size_t i = 0; i < m_interface.sidesA.size(); ++i)
204 assembler->setFixedDofs(m_interface.sidesB[i].patch,
205 m_interface.sidesB[i].side(),
206 disp.patch(m_interface.sidesA[i].patch).boundary(m_interface.sidesA[i].side())->coefs() -
207 ALEdisp.patch(m_interface.sidesB[i].patch).boundary(m_interface.sidesB[i].side())->coefs());
208 solverNL->reset();
209 solverNL->solve();
210 assembler->constructSolution(solverNL->solution(),solverNL->allFixedDofs(),ALEdisp);
211 if (m_options.getSwitch("Check"))
212 return checkDisplacement(assembler->patches(),ALEdisp);
213 else
214 return -1;
215}
216
217template <class T>
218void gsALE<T>::saveState()
219{
220 if (methodALE == ale_method::TINE || methodALE == ale_method::TINE_StVK)
221 solverNL->saveState();
222 ALEdispSaved.clear();
223 for (size_t p = 0; p < ALEdisp.nPatches(); ++p)
224 ALEdispSaved.addPatch(ALEdisp.patch(p).clone());
225 hasSavedState = true;
226}
227
228template <class T>
229void gsALE<T>::recoverState()
230{
231 GISMO_ENSURE(hasSavedState,"No state saved!");
232 if (methodALE == ale_method::TINE || methodALE == ale_method::TINE_StVK)
233 solverNL->recoverState();
234 if (methodALE == ale_method::IHE || methodALE == ale_method::ILE || methodALE == ale_method::IBHE)
235 for (size_t p = 0; p < ALEdisp.nPatches(); ++p)
236 assembler->patches().patch(p).coefs() += ALEdispSaved.patch(p).coefs() - ALEdisp.patch(p).coefs();
237 for (size_t p = 0; p < ALEdisp.nPatches(); ++p)
238 ALEdisp.patch(p).coefs() = ALEdispSaved.patch(p).coefs();
239}
240
241} // namespace ends
gsALE.h
Implementation of mesh deformation method for partitioned FSI solver.
index_t
#define index_t
Definition gsConfig.h:32
gsConstantFunction.h
Provides declaration of ConstantFunction class.
GISMO_ENSURE
#define GISMO_ENSURE(cond, message)
Definition gsDebug.h:102
gsElPoissonAssembler.h
Provides stiffness matrix for Poisson's equations.
gsElasticityAssembler.h
Provides linear and nonlinear elasticity systems for 2D plain strain and 3D continua.
gsGeoUtils.h
Provides isogeometric meshing and modelling routines.
gsIterative.h
A class providing an iterative solver for nonlinear problems.
gismo
The G+Smo namespace, containing all definitions for the library.
gismo::checkDisplacement
index_t checkDisplacement(gsMultiPatch< T > const &domain, gsMultiPatch< T > const &displacement)
Checks whether the deformed configuration is bijective, i.e. det(Jac(geo+disp)) > 0; returns -1 if ye...
Definition gsGeoUtils.hpp:245
gismo::checkGeometry
index_t checkGeometry(gsMultiPatch< T > const &domain)
Checks whether configuration is bijective, i.e. det(Jac(geo)) > 0; returns -1 if yes or the number of...
Definition gsGeoUtils.hpp:200
gsBiharmonicAssembler.h
Provides stiffness matrix for bi-harmonic equation in the mixed formulation.
gismo::ale_method::method
method
Definition gsBaseUtils.h:29