G+Smo  25.01.0
Geometry + Simulation Modules
 
Loading...
Searching...
No Matches
gsThinShellAssembler.hpp
Go to the documentation of this file.
1
16#pragma once
17
18#include <gsKLShell/src/gsThinShellAssembler.h>
19#include <gsKLShell/src/gsMaterialMatrixBase.h>
20#include <gsKLShell/src/gsMaterialMatrixIntegrate.h>
21#include <gsKLShell/src/gsMaterialMatrixEval.h>
22
23#include <gsPde/gsBoundaryConditions.h>
24
25#include <gsCore/gsFunctionExpr.h>
26#include <gsCore/gsConstantFunction.h>
27#include <gsCore/gsPiecewiseFunction.h>
28
29#include <unordered_set>
30
31namespace gismo
32{
33
34template<short_t d, class T, bool bending>
35gsThinShellAssembler<d, T, bending>::gsThinShellAssembler(const gsMultiPatch<T> & patches,
36 const gsMultiBasis<T> & basis,
37 const gsBoundaryConditions<T> & bconditions,
38 const gsFunctionSet<T> & surface_force,
39 const gsMaterialMatrixContainer<T> & materialMatrices
40 )
41 :
42 m_patches(patches),
43 m_basis(basis),
44 m_spaceBasis(&m_basis),
45 m_bcs(bconditions),
46 m_forceFun(&surface_force),
47 m_materialMatrices(materialMatrices)
48{
49 // surface forces defined in the parametric domain (ONLY WORKS FOR 3D)
50 m_parametricForce = (surface_force.domainDim()==2 && (d==2||d==3));
51
52 this->_defaultOptions();
53 this->_getOptions();
54 this->_initialize();
55}
56
57template<short_t d, class T, bool bending>
58gsThinShellAssembler<d, T, bending>::gsThinShellAssembler(const gsMultiPatch<T> & patches,
59 const gsMultiBasis<T> & basis,
60 const gsBoundaryConditions<T> & bconditions,
61 const gsFunctionSet<T> & surface_force,
62 typename gsMaterialMatrixBase<T>::uPtr & materialMatrix
63 )
64:
65gsThinShellAssembler<d, T, bending>(patches,basis,bconditions,surface_force,materialMatrix.get())
66{
67
68}
69
70template<short_t d, class T, bool bending>
71gsThinShellAssembler<d, T, bending>::gsThinShellAssembler(const gsMultiPatch<T> & patches,
72 const gsMultiBasis<T> & basis,
73 const gsBoundaryConditions<T> & bconditions,
74 const gsFunctionSet<T> & surface_force,
75 gsMaterialMatrixBase<T> * materialMatrix
76 )
77 :
78 m_patches(patches),
79 m_basis(basis),
80 m_spaceBasis(&basis),
81 m_bcs(bconditions),
82 m_forceFun(&surface_force)
83{
84 m_materialMatrices = gsMaterialMatrixContainer<T>(m_patches.nPatches());
85 GISMO_ASSERT(materialMatrix!=nullptr,"Material matrix is incomplete!");
86 GISMO_ASSERT(materialMatrix->initialized(),"Material matrix is incomplete!");
87 for (size_t p=0; p!=m_patches.nPatches(); p++)
88 m_materialMatrices.set(p,materialMatrix);
89
90 // surface forces defined in the parametric domain (ONLY WORKS FOR 3D)
91 m_parametricForce = (surface_force.domainDim()==2 && (d==2||d==3));
92
93 this->_defaultOptions();
94 this->_getOptions();
95 this->_initialize();
96}
97
98template<short_t d, class T, bool bending>
99gsThinShellAssembler<d, T, bending>& gsThinShellAssembler<d, T, bending>::operator=( const gsThinShellAssembler& other )
100{
101 if (this!=&other)
102 {
103 m_mapper=other.m_mapper;
104 m_assembler=other.m_assembler;
105 m_evaluator=other.m_evaluator;
106 m_patches=other.m_patches;
107 m_itpatches=other.m_itpatches;
108 m_basis=other.m_basis;
109 m_spaceBasis=other.m_spaceBasis;
110 m_bcs=other.m_bcs;
111 m_ddofs=other.m_ddofs;
112 m_mass=other.m_mass;
113 m_parametricForce=other.m_parametricForce;
114 m_forceFun=other.m_forceFun;
115 m_foundFun=other.m_foundFun;
116 m_pressFun=other.m_pressFun;
117 m_materialMatrices=other.m_materialMatrices;
118 m_pLoads=other.m_pLoads;
119 m_pMass=other.m_pMass;
120 m_solvector=other.m_solvector;
121 m_rhs=other.m_rhs;
122 m_options=other.m_options;
123 m_foundInd=other.m_foundInd;
124 m_pressInd=other.m_pressInd;
125 m_continuity=other.m_continuity;
126 m_alpha_d_bc=other.m_alpha_d_bc;
127 m_alpha_r_bc=other.m_alpha_r_bc;
128 m_alpha_d_ifc=other.m_alpha_d_ifc;
129 m_alpha_r_ifc=other.m_alpha_r_ifc;
130 m_IfcDefault=other.m_IfcDefault;
131 m_inPlane=other.m_inPlane;
132 m_outPlane=other.m_outPlane;
133 m_uncoupled=other.m_uncoupled;
134 m_strongC0=other.m_strongC0;
135 m_weakC0=other.m_weakC0;
136 m_strongC1=other.m_strongC1;
137 m_weakC1=other.m_weakC1;
138 m_unassigned=other.m_unassigned;
139
140 // To do: make copy constructor for the gsExprAssembler
141 m_assembler.setIntegrationElements(m_basis);
142 GISMO_ENSURE(m_options.hasGroup("ExprAssembler"),"The option list does not contain options with the label 'ExprAssembler'!");
143 m_assembler.setOptions(m_options.getGroup("ExprAssembler"));
144 }
145 return *this;
146}
147
148template<short_t d, class T, bool bending>
149gsThinShellAssembler<d, T, bending>& gsThinShellAssembler<d, T, bending>::operator=( gsThinShellAssembler&& other )
150{
151 m_mapper=give(other.m_mapper);
152 m_assembler=give(other.m_assembler);
153 m_evaluator=give(other.m_evaluator);
154 m_patches=give(other.m_patches);
155 m_itpatches=give(other.m_itpatches);
156 m_basis=give(other.m_basis);
157 m_spaceBasis=give(other.m_spaceBasis);
158 m_bcs=give(other.m_bcs);
159 m_ddofs=give(other.m_ddofs);
160 m_mass=give(other.m_mass);
161 m_parametricForce=give(other.m_parametricForce);
162 m_forceFun=give(other.m_forceFun);
163 m_foundFun=give(other.m_foundFun);
164 m_pressFun=give(other.m_pressFun);
165 m_materialMatrices=give(other.m_materialMatrices);
166 m_pLoads=give(other.m_pLoads);
167 m_pMass=give(other.m_pMass);
168 m_solvector=give(other.m_solvector);
169 m_rhs=give(other.m_rhs);
170 m_options=give(other.m_options);
171 m_foundInd=give(other.m_foundInd);
172 m_pressInd=give(other.m_pressInd);
173 m_continuity=give(other.m_continuity);
174 m_alpha_d_bc=give(other.m_alpha_d_bc);
175 m_alpha_r_bc=give(other.m_alpha_r_bc);
176 m_alpha_d_ifc=give(other.m_alpha_d_ifc);
177 m_alpha_r_ifc=give(other.m_alpha_r_ifc);
178 m_IfcDefault=give(other.m_IfcDefault);
179 m_inPlane=give(other.m_inPlane);
180 m_outPlane=give(other.m_outPlane);
181 m_uncoupled=give(other.m_uncoupled);
182 m_strongC0=give(other.m_strongC0);
183 m_weakC0=give(other.m_weakC0);
184 m_strongC1=give(other.m_strongC1);
185 m_weakC1=give(other.m_weakC1);
186 m_unassigned=give(other.m_unassigned);
187 // To do: make copy constructor for the gsExprAssembler
188 m_assembler.setIntegrationElements(m_basis);
189 GISMO_ENSURE(m_options.hasGroup("ExprAssembler"),"The option list does not contain options with the label 'ExprAssembler'!");
190 m_assembler.setOptions(m_options.getGroup("ExprAssembler"));
191 return *this;
192}
193
194template<short_t d, class T, bool bending>
195void gsThinShellAssembler<d, T, bending>::_defaultOptions()
196{
197 m_options.addReal("WeakDirichlet","Penalty parameter weak dirichlet conditions",1e3);
198 m_options.addReal("WeakClamped","Penalty parameter weak clamped conditions",1e3);
199 m_options.addInt("Continuity","Set the continuity for the space",-1);
200
201 m_options.addReal("IfcPenalty","Penalty parameter weak coupling conditions on the interface",1e3);
202 m_options.addInt("IfcDefault","Default weak(!) interface coupling; C^k, k={-1,0,1}",1);
203 m_options.addString("Solver","Sparse linear solver", "CGDiagonal");
204
205 // Assembler options
206 gsOptionList assemblerOptions = m_assembler.defaultOptions().wrapIntoGroup("ExprAssembler");
207 m_options.update(assemblerOptions,gsOptionList::addIfUnknown);
208
209 m_continuity = -1;
210}
211
212
213template <short_t d, class T, bool bending>
214void gsThinShellAssembler<d, T, bending>::_getOptions()
215{
216 // If the continuity changed, we need to re-initialize the space.
217 index_t continuity = m_continuity;
218 m_continuity = m_options.getInt("Continuity");
219 if (continuity != m_options.getInt("Continuity"))
220 this->_initialize();
221
222 m_alpha_d_bc = m_options.getReal("WeakDirichlet");
223 m_alpha_r_bc = m_options.getReal("WeakClamped");
224 m_alpha_d_ifc = m_alpha_r_ifc = m_options.getReal("IfcPenalty");
225 m_IfcDefault = m_options.getInt("IfcDefault");
226}
227
228template<short_t d, class T, bool bending>
229void gsThinShellAssembler<d, T, bending>::setOptions(gsOptionList & options)
230{
231 // Check if the continuity option changed
232 // Get old continuity
233 index_t continuity = m_options.getInt("Continuity");
234
235 m_options.update(options,gsOptionList::ignoreIfUnknown);
236
237 // If the continuity changed, we need to re-initialize the space.
238 if (continuity != m_options.getInt("Continuity"))
239 this->_initialize();
240}
241
242template<short_t d, class T, bool bending>
243void gsThinShellAssembler<d, T, bending>::_initialize()
244{
245 //gsInfo<<"Active options:\n"<< m_assembler.options() <<"\n";
246
247 // Elements used for numerical integration
248 m_assembler.setIntegrationElements(m_basis);
249 GISMO_ENSURE(m_options.hasGroup("ExprAssembler"),"The option list does not contain options with the label 'ExprAssembler'!");
250 m_assembler.setOptions(m_options.getGroup("ExprAssembler"));
251
252 GISMO_ASSERT(m_bcs.hasGeoMap(),"No geometry map was assigned to the boundary conditions. Use bc.setGeoMap to assign one!");
253
254 // Initialize the geometry maps
255 // geometryMap m_ori = m_assembler.getMap(m_patches);
256 // geometryMap m_def = m_assembler.getMap(*m_defpatches);
257
258 // Set the discretization space
259 space m_space = m_assembler.getSpace(*m_spaceBasis, d, 0); // last argument is the space ID
260
261 this->_assembleDirichlet();
262
263 m_ddofs = m_space.fixedPart();
264 m_mapper = m_space.mapper();
265
266 // foundation is off by default
267 m_foundInd = false;
268 // pressure is off by default
269 m_pressInd = false;
270
271 GISMO_ASSERT(m_forceFun->targetDim()==d,"Force must have " << d<<" dimensions but has "<<m_forceFun->targetDim());
272
273 // test interfaces on in-plane and out-of-plane connection and put them in respective containers
274 // _ifcTest();
275 // match interfaces where needed
276 // todo
277 // Put the interfaces in the right container depending on the
278
279}
280
281template<short_t d, class T, bool bending>
282void gsThinShellAssembler<d, T, bending>::addStrongC0(const gsBoxTopology::ifContainer & interfaces)
283{
284 m_strongC0 = interfaces;
285}
286
287template<short_t d, class T, bool bending>
288void gsThinShellAssembler<d, T, bending>::addStrongC1(const gsBoxTopology::ifContainer & interfaces)
289{
290 m_strongC1 = interfaces;
291}
292
293template<short_t d, class T, bool bending>
294void gsThinShellAssembler<d, T, bending>::addWeakC0(const gsBoxTopology::ifContainer & interfaces)
295{
296 m_weakC0 = interfaces;
297}
298
299template<short_t d, class T, bool bending>
300void gsThinShellAssembler<d, T, bending>::addWeakC1(const gsBoxTopology::ifContainer & interfaces)
301{
302 m_weakC1 = interfaces;
303}
304
305template<short_t d, class T, bool bending>
306void gsThinShellAssembler<d, T, bending>::addUncoupled(const gsBoxTopology::ifContainer & interfaces)
307{
308 m_uncoupled = interfaces;
309}
310
311template<short_t d, class T, bool bending>
312void gsThinShellAssembler<d, T, bending>::initInterfaces()
313{
314 this->_getOptions();
315 // Find unassigned interfaces and add them to the right containers
316 for (gsBoxTopology::const_iiterator it = m_patches.topology().iBegin(); it!=m_patches.topology().iEnd(); it++)
317 {
318 if (
319 std::find(m_strongC0.begin(), m_strongC0.end(), *it) == m_strongC0.end() // m_strongC0 does not contain *it
320 && std::find(m_strongC1.begin(), m_strongC1.end(), *it) == m_strongC1.end() // m_strongC1 does not contain *it
321 && std::find(m_weakC0.begin(), m_weakC0.end(), *it) == m_weakC0.end() // m_weakC0 does not contain *it
322 && std::find(m_weakC1.begin(), m_weakC1.end(), *it) == m_weakC1.end() // m_weakC1 does not contain *it
323 && std::find(m_uncoupled.begin(), m_uncoupled.end(), *it) == m_uncoupled.end() // m_uncoupled does not contain *it
324 )
325 {
326 if (m_IfcDefault==-1)
327 continue;
328 else if (m_IfcDefault==0)
329 m_weakC0.push_back(*it);
330 else if (m_IfcDefault==1)
331 m_weakC1.push_back(*it);
332 else
333 GISMO_ERROR("Option unknown");
334 }
335 }
336
337 // Set strong C0 using the setup function.
338}
339
340template<short_t d, class T, bool bending>
341void gsThinShellAssembler<d, T, bending>::_assembleNeumann()
342{
343 _assembleNeumann_impl<d>();
344}
345
346template <short_t d, class T, bool bending>
347template <short_t _d>
348typename std::enable_if<(_d==3), void>::type
349gsThinShellAssembler<d, T, bending>::_assembleNeumann_impl()
350{
351 geometryMap m_ori = m_assembler.getMap(m_patches);
352
353 space m_space = m_assembler.trialSpace(0); // last argument is the space ID
354 auto g_N = m_assembler.getBdrFunction(m_ori);
355 m_assembler.assembleBdr(m_bcs.get("Neumann"),m_space * g_N * meas(m_ori));
356}
357
358template <short_t d, class T, bool bending>
359template <short_t _d>
360typename std::enable_if<!(_d==3), void>::type
361gsThinShellAssembler<d, T, bending>::_assembleNeumann_impl()
362{
363 geometryMap m_ori = m_assembler.getMap(m_patches);
364
365 space m_space = m_assembler.trialSpace(0); // last argument is the space ID
366 auto g_N = m_assembler.getBdrFunction(m_ori);
367 m_assembler.assembleBdr(m_bcs.get("Neumann"), m_space * g_N * meas(m_ori));
368}
369
370template<short_t d, class T, bool bending>
371template<bool _matrix>
372void gsThinShellAssembler<d, T, bending>::_assemblePressure(const gsFunction<T> & pressFun)
373{
374 this->_getOptions();
375 _assemblePressure_impl<d,_matrix>(pressFun);
376}
377
378template <short_t d, class T, bool bending>
379template <short_t _d, bool matrix>
380typename std::enable_if<(_d==3) && matrix, void>::type
381gsThinShellAssembler<d, T, bending>::_assemblePressure_impl(const gsFunction<T> &)
382{
383 // No matrix contribution for the linear case
384}
385
386template <short_t d, class T, bool bending>
387template <short_t _d, bool _matrix>
388typename std::enable_if<(_d==3) && !_matrix, void>::type
389gsThinShellAssembler<d, T, bending>::_assemblePressure_impl(const gsFunction<T> & pressFun)
390{
391 gsMultiPatch<T> & defpatches = m_patches;
392 geometryMap m_ori = m_assembler.getMap(m_patches);
393 geometryMap m_def = m_assembler.getMap(defpatches);
394
395 space m_space = m_assembler.trialSpace(0); // last argument is the space ID
396
397 auto m_pressure = m_assembler.getCoeff(pressFun, m_ori);
398 GISMO_ASSERT(pressFun.targetDim()==1,"Pressure function has dimension "<<pressFun.targetDim()<<", but expected 1");
399
400 m_assembler.assemble(
401 m_pressure.val() * m_space * usn(m_def) * meas(m_ori)
402 );
403}
404
405template <short_t d, class T, bool bending>
406template <short_t _d, bool _matrix>
407typename std::enable_if<!(_d==3), void>::type
408gsThinShellAssembler<d, T, bending>::_assemblePressure_impl(const gsFunction<T> &)
409{
410 // Since pressure works out-of-plane, this function has no effect
411}
412
413template<short_t d, class T, bool bending>
414template<bool _matrix>
415void gsThinShellAssembler<d, T, bending>::_assemblePressure(const gsFunction<T> & pressFun, const gsFunctionSet<T> & deformed)
416{
417 this->_getOptions();
418 _assemblePressure_impl<d,_matrix>(pressFun,deformed);
419}
420
421// assembles eq 3.26 from http://resolver.tudelft.nl/uuid:56c0cc91-643d-4817-9702-93fedce5fd78
422template <short_t d, class T, bool bending>
423template <short_t _d, bool matrix>
424typename std::enable_if<(_d==3) && matrix, void>::type
425gsThinShellAssembler<d, T, bending>::_assemblePressure_impl(const gsFunction<T> & pressFun, const gsFunctionSet<T> & deformed)
426{
427 geometryMap m_ori = m_assembler.getMap(m_patches);
428 geometryMap m_def = m_assembler.getMap(deformed);
429
430 space m_space = m_assembler.trialSpace(0); // last argument is the space ID
431
432 auto m_pressure = m_assembler.getCoeff(pressFun, m_ori);
433 GISMO_ASSERT(pressFun.targetDim()==1,"Pressure function has dimension "<<pressFun.targetDim()<<", but expected 1");
434
435 m_assembler.assemble(
436 -m_pressure.val() * m_space * var1(m_space,m_def).tr()* meas(m_ori)
437 //-m_pressure.val() * jac(m_space) * sn(m_def).normalized() * meas(m_ori)
438 );
439}
440
441// assembles eq 3.25 from http://resolver.tudelft.nl/uuid:56c0cc91-643d-4817-9702-93fedce5fd78
442template <short_t d, class T, bool bending>
443template <short_t _d, bool _matrix>
444typename std::enable_if<(_d==3) && !_matrix, void>::type
445gsThinShellAssembler<d, T, bending>::_assemblePressure_impl(const gsFunction<T> & pressFun, const gsFunctionSet<T> & deformed)
446{
447 geometryMap m_ori = m_assembler.getMap(m_patches);
448 geometryMap m_def = m_assembler.getMap(deformed);
449
450 space m_space = m_assembler.trialSpace(0); // last argument is the space ID
451
452 auto m_pressure = m_assembler.getCoeff(pressFun, m_ori);
453 GISMO_ASSERT(pressFun.targetDim()==1,"Pressure function has dimension "<<pressFun.targetDim()<<", but expected 1");
454
455 // Assemble vector
456 m_assembler.assemble(
457 m_pressure.val() * m_space * sn(m_def).normalized() * meas(m_ori)
458 );
459}
460
461template <short_t d, class T, bool bending>
462template <short_t _d, bool _matrix>
463typename std::enable_if<!(_d==3), void>::type
464gsThinShellAssembler<d, T, bending>::_assemblePressure_impl(const gsFunction<T> & /*presFun*/, const gsFunctionSet<T> & /*deformed*/)
465{
466 // Since pressure works out-of-plane, this function has no effect
467}
468
469template<short_t d, class T, bool bending>
470template<bool _matrix>
471void gsThinShellAssembler<d, T, bending>::_assembleFoundation(const gsFunction<T> & foundFun)
472{
473 this->_getOptions();
474 _assembleFoundation_impl<d,_matrix>(foundFun);
475}
476
477template <short_t d, class T, bool bending>
478template <short_t _d, bool _matrix>
479typename std::enable_if<(_d==3) && _matrix, void>::type
480gsThinShellAssembler<d, T, bending>::_assembleFoundation_impl(const gsFunction<T> & foundFun)
481{
482 geometryMap m_ori = m_assembler.getMap(m_patches);
483
484 space m_space = m_assembler.trialSpace(0); // last argument is the space ID
485
486 auto m_foundation = m_assembler.getCoeff(foundFun, m_ori);
487 GISMO_ASSERT(foundFun.targetDim()==3,"Foundation function has dimension "<<foundFun.targetDim()<<", but expected 3");
488
489 m_assembler.assemble(
490 m_space * m_foundation.asDiag() * m_space.tr() * meas(m_ori)
491 );
492}
493
494// assembles eq 3.27 from http://resolver.tudelft.nl/uuid:56c0cc91-643d-4817-9702-93fedce5fd78
495template <short_t d, class T, bool bending>
496template <short_t _d, bool _matrix>
497typename std::enable_if<(_d==3) && !_matrix, void>::type
498gsThinShellAssembler<d, T, bending>::_assembleFoundation_impl(const gsFunction<T> & /* foundFun */)
499{
500 // No rhs contribution for the linear case
501}
502
503template <short_t d, class T, bool bending>
504template <short_t _d, bool _matrix>
505typename std::enable_if<!(_d==3), void>::type
506gsThinShellAssembler<d, T, bending>::_assembleFoundation_impl(const gsFunction<T> & )
507{
508 // Since foundation works out-of-plane, this function has no effect
509}
510
511template<short_t d, class T, bool bending>
512template<bool _matrix>
513void gsThinShellAssembler<d, T, bending>::_assembleFoundation(const gsFunction<T> & foundFun, const gsFunctionSet<T> & deformed)
514{
515 this->_getOptions();
516 _assembleFoundation_impl<d,_matrix>(foundFun,deformed);
517}
518
519// assembles eq 3.28 from http://resolver.tudelft.nl/uuid:56c0cc91-643d-4817-9702-93fedce5fd78
520template <short_t d, class T, bool bending>
521template <short_t _d, bool matrix>
522typename std::enable_if<(_d==3) && matrix, void>::type
523gsThinShellAssembler<d, T, bending>::_assembleFoundation_impl(const gsFunction<T> & foundFun, const gsFunctionSet<T> & /*deformed*/)
524{
525 geometryMap m_ori = m_assembler.getMap(m_patches);
526
527 space m_space = m_assembler.trialSpace(0); // last argument is the space ID
528
529 auto m_foundation = m_assembler.getCoeff(foundFun, m_ori);
530 GISMO_ASSERT(foundFun.targetDim()==3,"Foundation function has dimension "<<foundFun.targetDim()<<", but expected 3");
531
532 m_assembler.assemble(
533 m_space * m_foundation.asDiag() * m_space.tr() * meas(m_ori)
534 );
535}
536
537template <short_t d, class T, bool bending>
538template <short_t _d, bool _matrix>
539typename std::enable_if<(_d==3) && !_matrix, void>::type
540gsThinShellAssembler<d, T, bending>::_assembleFoundation_impl(const gsFunction<T> & foundFun, const gsFunctionSet<T> & deformed)
541{
542 geometryMap m_ori = m_assembler.getMap(m_patches);
543 geometryMap m_def = m_assembler.getMap(deformed);
544
545 space m_space = m_assembler.trialSpace(0); // last argument is the space ID
546
547 auto m_foundation = m_assembler.getCoeff(foundFun, m_ori);
548 GISMO_ASSERT(foundFun.targetDim()==3,"Foundation function has dimension "<<foundFun.targetDim()<<", but expected 3");
549
550 // Assemble vector
551 m_assembler.assemble(
552 m_space * m_foundation.asDiag() * (m_def - m_ori) * meas(m_ori) // [v_x,v_y,v_z] diag([k_x,k_y,k_z]) [u_x; u_y; u_z]
553 );
554}
555
556template <short_t d, class T, bool bending>
557template <short_t _d, bool _matrix>
558typename std::enable_if<!(_d==3), void>::type
559gsThinShellAssembler<d, T, bending>::_assembleFoundation_impl(const gsFunction<T> & , const gsFunctionSet<T> & )
560{
561 // Since foundation works out-of-plane, this function has no effect
562}
563
564template<short_t d, class T, bool bending>
565template<bool _matrix>
566void gsThinShellAssembler<d, T, bending>::_assembleWeakBCs()
567{
568 this->_getOptions();
569 _assembleWeakBCs_impl<d,_matrix>();
570}
571
572template <short_t d, class T, bool bending>
573template <short_t _d, bool _matrix>
574typename std::enable_if<(_d==3) && _matrix, void>::type
575gsThinShellAssembler<d, T, bending>::_assembleWeakBCs_impl()
576{
577 gsMultiPatch<T> & defpatches = m_patches;
578 geometryMap m_ori = m_assembler.getMap(m_patches);
579
580 space m_space = m_assembler.trialSpace(0); // last argument is the space ID
581 // auto g_N = m_assembler.getBdrFunction(m_ori);
582
583 gsMaterialMatrixIntegrate<T,MaterialOutput::MatrixA> m_mmA(m_materialMatrices,&m_patches,&defpatches);
584 gsMaterialMatrixIntegrate<T,MaterialOutput::MatrixD> m_mmD(m_materialMatrices,&m_patches,&defpatches);
585 auto mmA = m_assembler.getCoeff(m_mmA);
586 auto mmD = m_assembler.getCoeff(m_mmD);
587
588 auto cart2cov = cartcov(m_ori);
589 auto con2cartI = cartcon(m_ori).inv();
590
591 auto mmAcart = (con2cartI * reshape(mmA,3,3) * cart2cov);
592 auto mmDcart = (con2cartI * reshape(mmD,3,3) * cart2cov);
593
594 element el = m_assembler.getElement();
595 auto alpha_d = m_alpha_d_bc * reshape(mmAcart,9,1).max().val() / el.area(m_ori);
596 auto alpha_r = m_alpha_r_bc * reshape(mmDcart,9,1).max().val() / el.area(m_ori);
597
598 // Weak BCs
599 m_assembler.assembleBdr
600 (
601 m_bcs.get("Weak Dirichlet")
602 ,
603 -(alpha_d * m_space * m_space.tr()) * meas(m_ori)
604 );
605
606 // for weak clamped
607 m_assembler.assembleBdr
608 (
609 m_bcs.get("Weak Clamped")
610 ,
611 (
612 alpha_r * ( ( var1(m_space,m_ori) * unv(m_ori) ) * ( var1(m_space,m_ori) * unv(m_ori) ).tr() )
613 ) * meas(m_ori)
614 );
615}
616
617template <short_t d, class T, bool bending>
618template <short_t _d, bool _matrix>
619typename std::enable_if<(_d==3) && !_matrix, void>::type
620gsThinShellAssembler<d, T, bending>::_assembleWeakBCs_impl()
621{
622 gsMultiPatch<T> & defpatches = m_patches;
623 geometryMap m_ori = m_assembler.getMap(m_patches);
624
625 space m_space = m_assembler.trialSpace(0); // last argument is the space ID
626 auto g_N = m_assembler.getBdrFunction(m_ori);
627
628 gsMaterialMatrixIntegrate<T,MaterialOutput::MatrixA> m_mmA(m_materialMatrices,&m_patches,&defpatches);
629 auto mmA = m_assembler.getCoeff(m_mmA);
630
631 auto cart2cov = cartcov(m_ori);
632 auto con2cartI = cartcon(m_ori).inv();
633
634 auto mmAcart = (con2cartI * reshape(mmA,3,3) * cart2cov);
635
636 element el = m_assembler.getElement();
637 auto alpha_d = m_alpha_d_bc * reshape(mmAcart,9,1).max().val() / el.area(m_ori);
638
639 // Weak BCs
640
641 m_assembler.assembleBdr
642 (
643 m_bcs.get("Weak Dirichlet")
644 ,
645 -(alpha_d * m_space * g_N ) * meas(m_ori)
646 );
647
648 // for weak clamped
649 // do nothing
650}
651
652template <short_t d, class T, bool bending>
653template <short_t _d, bool _matrix>
654typename std::enable_if<!(_d==3) && _matrix, void>::type
655gsThinShellAssembler<d, T, bending>::_assembleWeakBCs_impl()
656{
657 gsMultiPatch<T> & defpatches = m_patches;
658 geometryMap m_ori = m_assembler.getMap(m_patches);
659
660 space m_space = m_assembler.trialSpace(0); // last argument is the space ID
661 // auto g_N = m_assembler.getBdrFunction(m_ori);
662
663 gsMaterialMatrixIntegrate<T,MaterialOutput::MatrixA> m_mmA(m_materialMatrices,&m_patches,&defpatches);
664 auto mmA = m_assembler.getCoeff(m_mmA);
665
666 auto cart2cov = cartcov(m_ori);
667 auto con2cartI = cartcon(m_ori).inv();
668
669 auto mmAcart = (con2cartI * reshape(mmA,3,3) * cart2cov);
670
671 element el = m_assembler.getElement();
672 auto alpha_d = m_alpha_d_bc * reshape(mmAcart,9,1).max().val() / el.area(m_ori);
673
674 // Weak BCs
675 m_assembler.assembleBdr
676 (
677 m_bcs.get("Weak Dirichlet")
678 ,
679 -(alpha_d * m_space * m_space.tr()) * meas(m_ori)
680 );
681}
682
683template <short_t d, class T, bool bending>
684template <short_t _d, bool _matrix>
685typename std::enable_if<!(_d==3) && !_matrix, void>::type
686gsThinShellAssembler<d, T, bending>::_assembleWeakBCs_impl()
687{
688 geometryMap m_ori = m_assembler.getMap(m_patches);
689
690 space m_space = m_assembler.trialSpace(0); // last argument is the space ID
691 auto g_N = m_assembler.getBdrFunction(m_ori);
692
693 // Weak BCs
694 m_assembler.assembleBdr
695 (
696 m_bcs.get("Weak Dirichlet")
697 ,
698 (m_alpha_d_bc * (m_space * (m_ori - m_ori) - m_space * (g_N) )) * meas(m_ori)
699 );
700}
701
702template <short_t d, class T, bool bending>
703template <bool _matrix>
704void gsThinShellAssembler<d, T, bending>::_assembleWeakBCs(const gsFunctionSet<T> & deformed)
705{
706 this->_getOptions();
707 _assembleWeakBCs_impl<d,_matrix>(deformed);
708}
709
710template <short_t d, class T, bool bending>
711template <short_t _d, bool _matrix>
712typename std::enable_if<(_d==3) && _matrix, void>::type
713gsThinShellAssembler<d, T, bending>::_assembleWeakBCs_impl(const gsFunctionSet<T> & deformed)
714{
715 geometryMap m_ori = m_assembler.getMap(m_patches);
716 geometryMap m_def = m_assembler.getMap(deformed);
717
718 space m_space = m_assembler.trialSpace(0); // last argument is the space ID
719 // auto g_N = m_assembler.getBdrFunction(m_ori);
720
721 gsMaterialMatrixIntegrate<T,MaterialOutput::MatrixA> m_mmA(m_materialMatrices,&m_patches,&deformed);
722 gsMaterialMatrixIntegrate<T,MaterialOutput::MatrixD> m_mmD(m_materialMatrices,&m_patches,&deformed);
723 auto mmA = m_assembler.getCoeff(m_mmA);
724 auto mmD = m_assembler.getCoeff(m_mmD);
725
726 auto cart2cov = cartcov(m_ori);
727 auto con2cartI = cartcon(m_ori).inv();
728
729 auto mmAcart = (con2cartI * reshape(mmA,3,3) * cart2cov);
730 auto mmDcart = (con2cartI * reshape(mmD,3,3) * cart2cov);
731
732 element el = m_assembler.getElement();
733 auto alpha_d = m_alpha_d_bc * reshape(mmAcart,9,1).max().val() / el.area(m_ori);
734 auto alpha_r = m_alpha_r_bc * reshape(mmDcart,9,1).max().val() / el.area(m_ori);
735
736
737 auto du = m_def - m_ori;
738 auto dnN = ( usn(m_def).tr()*unv(m_ori) - usn(m_ori).tr()*unv(m_ori) ).val();
739
740 // Weak BCs
741 m_assembler.assembleBdr
742 (
743 m_bcs.get("Weak Dirichlet")
744 ,
745 -alpha_d * m_space * m_space.tr() * meas(m_ori)
746 );
747
748 // for weak clamped
749 m_assembler.assembleBdr
750 (
751 m_bcs.get("Weak Clamped")
752 ,
753 (
754 alpha_r * dnN * ( var2deriv2(m_space,m_space,m_def,unv(m_ori).tr()) )
755 +
756 alpha_r * ( ( var1(m_space,m_def) * unv(m_ori) ) * ( var1(m_space,m_def) * unv(m_ori) ).tr() )
757 ) * meas(m_ori)
758 );
759}
760
761template <short_t d, class T, bool bending>
762template <short_t _d, bool _matrix>
763typename std::enable_if<(_d==3) && !_matrix, void>::type
764gsThinShellAssembler<d, T, bending>::_assembleWeakBCs_impl(const gsFunctionSet<T> & deformed)
765{
766 geometryMap m_ori = m_assembler.getMap(m_patches);
767 geometryMap m_def = m_assembler.getMap(deformed);
768
769 space m_space = m_assembler.trialSpace(0); // last argument is the space ID
770 auto g_N = m_assembler.getBdrFunction(m_ori);
771
772 gsMaterialMatrixIntegrate<T,MaterialOutput::MatrixA> m_mmA(m_materialMatrices,&m_patches,&deformed);
773 gsMaterialMatrixIntegrate<T,MaterialOutput::MatrixD> m_mmD(m_materialMatrices,&m_patches,&deformed);
774 auto mmA = m_assembler.getCoeff(m_mmA);
775 auto mmD = m_assembler.getCoeff(m_mmD);
776
777 auto cart2cov = cartcov(m_ori);
778 auto con2cartI = cartcon(m_ori).inv();
779
780 auto mmAcart = (con2cartI * reshape(mmA,3,3) * cart2cov);
781 auto mmDcart = (con2cartI * reshape(mmD,3,3) * cart2cov);
782
783 element el = m_assembler.getElement();
784 auto alpha_d = m_alpha_d_bc * reshape(mmAcart,9,1).max().val() / el.area(m_ori);
785 auto alpha_r = m_alpha_r_bc * reshape(mmDcart,9,1).max().val() / el.area(m_ori);
786
787 auto du = m_def - m_ori;
788 auto dnN = ( usn(m_def).tr()*nv(m_ori) - usn(m_ori).tr()*nv(m_ori) ).val();
789
790 // Weak BCs
791 m_assembler.assembleBdr
792 (
793 m_bcs.get("Weak Dirichlet")
794 ,
795 alpha_d * (m_space * du - m_space * (g_N) ) * meas(m_ori)
796 );
797
798 // for weak clamped
799 m_assembler.assembleBdr
800 (
801 m_bcs.get("Weak Clamped")
802 ,
803 (
804 - alpha_r * dnN * ( var1(m_space,m_def) * usn(m_ori) )
805 ) * meas(m_ori)
806 );
807}
808
809template <short_t d, class T, bool bending>
810template <short_t _d, bool _matrix>
811typename std::enable_if<!(_d==3) && _matrix, void>::type
812gsThinShellAssembler<d, T, bending>::_assembleWeakBCs_impl(const gsFunctionSet<T> & deformed)
813{
814 geometryMap m_ori = m_assembler.getMap(m_patches);
815
816 space m_space = m_assembler.trialSpace(0); // last argument is the space ID
817
818 gsMaterialMatrixIntegrate<T,MaterialOutput::MatrixA> m_mmA(m_materialMatrices,&m_patches,&deformed);
819 auto mmA = m_assembler.getCoeff(m_mmA);
820
821 auto cart2cov = cartcov(m_ori);
822 auto con2cartI = cartcon(m_ori).inv();
823
824 auto mmAcart = (con2cartI * reshape(mmA,3,3) * cart2cov);
825
826 element el = m_assembler.getElement();
827 auto alpha_d = m_alpha_d_bc * reshape(mmAcart,9,1).max().val() / el.area(m_ori);
828
829 // Weak BCs
830 m_assembler.assembleBdr
831 (
832 m_bcs.get("Weak Dirichlet")
833 ,
834 -alpha_d * m_space * m_space.tr() * meas(m_ori)
835 );
836}
837
838template <short_t d, class T, bool bending>
839template <short_t _d, bool _matrix>
840typename std::enable_if<!(_d==3) && !_matrix, void>::type
841gsThinShellAssembler<d, T, bending>::_assembleWeakBCs_impl(const gsFunctionSet<T> & deformed)
842{
843 geometryMap m_ori = m_assembler.getMap(m_patches);
844 geometryMap m_def = m_assembler.getMap(deformed);
845
846 space m_space = m_assembler.trialSpace(0); // last argument is the space ID
847 auto g_N = m_assembler.getBdrFunction(m_ori);
848
849 gsMaterialMatrixIntegrate<T,MaterialOutput::MatrixA> m_mmA(m_materialMatrices,&m_patches,&deformed);
850 auto mmA = m_assembler.getCoeff(m_mmA);
851
852 auto cart2cov = cartcov(m_ori);
853 auto con2cartI = cartcon(m_ori).inv();
854
855 auto mmAcart = (con2cartI * reshape(mmA,3,3) * cart2cov);
856
857 element el = m_assembler.getElement();
858 auto alpha_d = m_alpha_d_bc * reshape(mmAcart,9,1).max().val() / el.area(m_ori);
859
860 // Weak BCs
861 m_assembler.assembleBdr
862 (
863 m_bcs.get("Weak Dirichlet")
864 ,
865 alpha_d * (m_space * (m_def - m_ori) - m_space * (g_N) ) * meas(m_ori)
866 );
867}
868
869template<short_t d, class T, bool bending>
870template<bool _matrix>
871void gsThinShellAssembler<d, T, bending>::_assembleWeakIfc()
872{
873 this->_getOptions();
874 if (m_weakC0.size()==0 && m_weakC1.size()==0)
875 return;
876 _assembleWeakIfc_impl<d,_matrix>();
877}
878
879template <short_t d, class T, bool bending>
880template <short_t _d, bool _matrix>
881typename std::enable_if<(_d==3) && _matrix, void>::type
882gsThinShellAssembler<d, T, bending>::_assembleWeakIfc_impl()
883{
884 gsMultiPatch<T> & defpatches = m_patches;
885 geometryMap m_ori = m_assembler.getMap(m_patches);
886
887 space m_space = m_assembler.trialSpace(0); // last argument is the space ID
888 // auto g_N = m_assembler.getBdrFunction(m_ori);
889
890 gsMaterialMatrixIntegrate<T,MaterialOutput::MatrixA> m_mmA(m_materialMatrices,&m_patches,&defpatches);
891 gsMaterialMatrixIntegrate<T,MaterialOutput::MatrixD> m_mmD(m_materialMatrices,&m_patches,&defpatches);
892 auto mmA = m_assembler.getCoeff(m_mmA);
893 auto mmD = m_assembler.getCoeff(m_mmD);
894
895 auto cart2cov = cartcov(m_ori);
896 auto con2cartI = cartcon(m_ori).inv();
897
898 auto mmAcart = (con2cartI * reshape(mmA,3,3) * cart2cov);
899 auto mmDcart = (con2cartI * reshape(mmD,3,3) * cart2cov);
900
901 element el = m_assembler.getElement();
902 auto h = (el.area(m_ori.left()) + el.area(m_ori.right())) / 2;
903 auto alpha_d = m_alpha_d_ifc * reshape(mmAcart,9,1).max().val() / h;
904 auto alpha_r = m_alpha_r_ifc * reshape(mmDcart,9,1).max().val() / h;
905
906 // C^0 coupling
907 m_assembler.assembleIfc(m_weakC0,
908 alpha_d * m_space.left() * m_space.left().tr() * meas(m_ori)
909 ,
910 -alpha_d * m_space.right()* m_space.left() .tr() * meas(m_ori)
911 ,
912 -alpha_d * m_space.left() * m_space.right().tr() * meas(m_ori)
913 ,
914 alpha_d * m_space.right()* m_space.right().tr() * meas(m_ori)
915 );
916
917 // C^1 coupling
918 // Penalty of out-of-plane coupling
919 // dW^pr / du_r --> second line
920 m_assembler.assembleIfc(m_weakC1,
921 alpha_r * ( var1(m_space.left(),m_ori.left()) * usn(m_ori.right()) ) * ( var1(m_space.left(),m_ori.left()) * usn(m_ori.right()) ).tr() * meas(m_ori) // left left
922 ,
923 alpha_r * ( var1(m_space.left(),m_ori.left()) * usn(m_ori.right()) ) * ( var1(m_space.right(),m_ori.right()) * usn(m_ori.left()) ).tr() * meas(m_ori) // left right
924 ,
925 alpha_r * ( var1(m_space.right(),m_ori.right()) * usn(m_ori.left()) ) * ( var1(m_space.left(),m_ori.left()) * usn(m_ori.right()) ).tr() * meas(m_ori) // right left
926 ,
927 alpha_r * ( var1(m_space.right(),m_ori.right()) * usn(m_ori.left()) ) * ( var1(m_space.right(),m_ori.right()) * usn(m_ori.left()) ).tr() * meas(m_ori) // right right
928 ,
929 // Symmetry
930 alpha_r * ( var1(m_space.right(),m_ori.right()) * usn(m_ori.left()) ) * ( var1(m_space.right(),m_ori.right()) * usn(m_ori.left()) ).tr() * meas(m_ori) // right right
931 ,
932 alpha_r * ( var1(m_space.right(),m_ori.right()) * usn(m_ori.left()) ) * ( var1(m_space.left(),m_ori.left()) * usn(m_ori.right()) ).tr() * meas(m_ori) // right left
933 ,
934 alpha_r * ( var1(m_space.left(),m_ori.left()) * usn(m_ori.right()) ) * ( var1(m_space.right(),m_ori.right()) * usn(m_ori.left()) ).tr() * meas(m_ori) // left right
935 ,
936 alpha_r * ( var1(m_space.left(),m_ori.left()) * usn(m_ori.right()) ) * ( var1(m_space.left(),m_ori.left()) * usn(m_ori.right()) ).tr() * meas(m_ori) // left left
937 );
938
939 // Penalty of in-plane coupling
940 // dW^pr / du_r --> fourth line
941 m_assembler.assembleIfc(m_weakC1,
942 alpha_r * ( ovar1(m_space.left() ,m_ori.left() ) * usn(m_ori.right()) ) * ( ovar1(m_space.left() ,m_ori.left() ) * usn(m_ori.right()) ).tr() * meas(m_ori) // left left
943 + // Symmetry
944 alpha_r * ( var1(m_space.left() ,m_ori.left() ) * unv(m_ori.right()) ) * ( var1(m_space.left() ,m_ori.left() ) * unv(m_ori.right()) ).tr() * meas(m_ori) // left left
945 ,
946 alpha_r * ( ovar1(m_space.left() ,m_ori.left() ) * usn(m_ori.right()) ) * ( var1(m_space.right(),m_ori.right()) * unv(m_ori.left() ) ).tr() * meas(m_ori) // left right
947 + // Symmetry
948 alpha_r * ( var1(m_space.left() ,m_ori.left() ) * unv(m_ori.right()) ) * ( ovar1(m_space.right(),m_ori.right()) * usn(m_ori.left() ) ).tr() * meas(m_ori) // left right
949 ,
950 alpha_r * ( var1(m_space.right(),m_ori.right()) * unv(m_ori.left() ) ) * ( ovar1(m_space.left() ,m_ori.left() ) * usn(m_ori.right()) ).tr() * meas(m_ori) // right left
951 + // Symmetry
952 alpha_r * ( ovar1(m_space.right(),m_ori.right()) * usn(m_ori.left() ) ) * ( var1(m_space.left() ,m_ori.left() ) * unv(m_ori.right()) ).tr() * meas(m_ori) // right left
953 ,
954 alpha_r * ( var1(m_space.right(),m_ori.right()) * unv(m_ori.left() ) ) * ( var1(m_space.right(),m_ori.right()) * unv(m_ori.left() ) ).tr() * meas(m_ori) // right right
955 + // Symmetry
956 alpha_r * ( ovar1(m_space.right(),m_ori.right()) * usn(m_ori.left() ) ) * ( ovar1(m_space.right(),m_ori.right()) * usn(m_ori.left() ) ).tr() * meas(m_ori) // right right
957 );
958}
959
960template <short_t d, class T, bool bending>
961template <short_t _d, bool _matrix>
962typename std::enable_if<(_d==3) && !_matrix, void>::type
963gsThinShellAssembler<d, T, bending>::_assembleWeakIfc_impl()
964{
965/*
966 empty
967 */
968}
969
970template <short_t d, class T, bool bending>
971template <short_t _d, bool _matrix>
972typename std::enable_if<!(_d==3) && _matrix, void>::type
973gsThinShellAssembler<d, T, bending>::_assembleWeakIfc_impl()
974{
975 gsMultiPatch<T> & defpatches = m_patches;
976 geometryMap m_ori = m_assembler.getMap(m_patches);
977
978 space m_space = m_assembler.trialSpace(0); // last argument is the space ID
979
980 gsMaterialMatrixIntegrate<T,MaterialOutput::MatrixA> m_mmA(m_materialMatrices,&m_patches,&defpatches);
981 auto mmA = m_assembler.getCoeff(m_mmA);
982
983 auto cart2cov = cartcov(m_ori);
984 auto con2cartI = cartcon(m_ori).inv();
985
986 auto mmAcart = (con2cartI * reshape(mmA,3,3) * cart2cov);
987
988 element el = m_assembler.getElement();
989 auto h = (el.area(m_ori.left()) + el.area(m_ori.right())) / 2;
990 auto alpha_d = m_alpha_d_ifc * reshape(mmAcart,9,1).max().val() / h;
991
992 // C^0 coupling
993 m_assembler.assembleIfc(m_weakC0,
994 alpha_d * m_space.left() * m_space.left().tr() * meas(m_ori) * meas(m_ori)
995 ,
996 -alpha_d * m_space.right()* m_space.left() .tr() * meas(m_ori) * meas(m_ori)
997 ,
998 -alpha_d * m_space.left() * m_space.right().tr() * meas(m_ori) * meas(m_ori)
999 ,
1000 alpha_d * m_space.right()* m_space.right().tr() * meas(m_ori) * meas(m_ori)
1001 );
1002
1003 // C^1 coupling DOES NOT CONTRIBUTE IN 2D PROBLEMS
1004}
1005
1006template <short_t d, class T, bool bending>
1007template <short_t _d, bool _matrix>
1008typename std::enable_if<!(_d==3) && !_matrix, void>::type
1009gsThinShellAssembler<d, T, bending>::_assembleWeakIfc_impl()
1010{
1011/*
1012 empty
1013 */
1014}
1015
1016template<short_t d, class T, bool bending>
1017template<bool _matrix>
1018void gsThinShellAssembler<d, T, bending>::_assembleWeakIfc(const gsFunctionSet<T> & deformed)
1019{
1020 this->_getOptions();
1021 _assembleWeakIfc_impl<d,_matrix>(deformed);
1022}
1023
1024template <short_t d, class T, bool bending>
1025template <short_t _d, bool _matrix>
1026typename std::enable_if<(_d==3) && _matrix, void>::type
1027gsThinShellAssembler<d, T, bending>::_assembleWeakIfc_impl(const gsFunctionSet<T> & deformed)
1028{
1029 geometryMap m_ori = m_assembler.getMap(m_patches);
1030 geometryMap m_def = m_assembler.getMap(deformed);
1031
1032 space m_space = m_assembler.trialSpace(0); // last argument is the space ID
1033
1034 gsMaterialMatrixIntegrate<T,MaterialOutput::MatrixA> m_mmA(m_materialMatrices,&m_patches,&deformed);
1035 gsMaterialMatrixIntegrate<T,MaterialOutput::MatrixD> m_mmD(m_materialMatrices,&m_patches,&deformed);
1036 auto mmA = m_assembler.getCoeff(m_mmA);
1037 auto mmD = m_assembler.getCoeff(m_mmD);
1038
1039 auto cart2cov = cartcov(m_ori);
1040 auto con2cartI = cartcon(m_ori).inv();
1041
1042 auto mmAcart = (con2cartI * reshape(mmA,3,3) * cart2cov);
1043 auto mmDcart = (con2cartI * reshape(mmD,3,3) * cart2cov);
1044
1045 element el = m_assembler.getElement();
1046 auto h = (el.area(m_ori.left()) + el.area(m_ori.right())) / 2;
1047 auto alpha_d = m_alpha_d_ifc * reshape(mmAcart,9,1).max().val() / h;
1048 auto alpha_r = m_alpha_r_ifc * reshape(mmDcart,9,1).max().val() / h;
1049
1050 auto du = ((m_def.left()-m_ori.left()) - (m_def.right()-m_ori.right()));
1051
1052 auto dN_lr = (usn(m_def.left()).tr()*usn(m_def.right())
1053 - usn(m_ori.left()).tr()*usn(m_ori.right())).val();
1054
1055 auto dN_rl = (usn(m_def.right()).tr()*usn(m_def.left())
1056 - usn(m_ori.right()).tr()*usn(m_ori.left())).val();
1057
1058 auto dnN_lr= (unv(m_def.left()).tr()*usn(m_def.right())
1059 - unv(m_ori.left()).tr()*usn(m_ori.right())).val();
1060
1061 auto dnN_rl= (unv(m_def.right()).tr()*usn(m_def.left())
1062 - unv(m_ori.right()).tr()*usn(m_ori.left())).val();
1063
1064 // C^0 coupling
1065 m_assembler.assembleIfc(m_weakC0,
1066 alpha_d * m_space.left() * m_space.left().tr() * meas(m_ori)
1067 ,
1068 -alpha_d * m_space.right()* m_space.left() .tr() * meas(m_ori)
1069 ,
1070 -alpha_d * m_space.left() * m_space.right().tr() * meas(m_ori)
1071 ,
1072 alpha_d * m_space.right()* m_space.right().tr() * meas(m_ori)
1073 );
1074
1075 // C^1 coupling
1076 // Penalty of out-of-plane coupling
1077 // dW^pr / du_r --> first line
1078 m_assembler.assembleIfc(m_weakC1,
1079 alpha_r * dN_lr * var2(m_space.left() ,m_space.left() ,m_def.left() ,usn(m_def.right()).tr() ) * meas(m_ori) // left left
1080 +//Symmetry
1081 alpha_r * dN_rl * var2( m_space.left(),m_space.left(),m_def.left(),usn(m_def.right() ).tr() ) * meas(m_ori) // left left
1082 ,
1083 alpha_r * dN_lr * ( var1(m_space.left() ,m_def.left() ) * var1(m_space.right(),m_def.right()).tr() ) * meas(m_ori)// left right
1084 +//Symmetry
1085 alpha_r * dN_rl * ( var1(m_space.left(),m_def.left()) * var1(m_space.right() ,m_def.right() ).tr() ) * meas(m_ori)// left right
1086 ,
1087 alpha_r * dN_lr * ( var1(m_space.right(),m_def.right()) * var1(m_space.left() ,m_def.left() ).tr() ) * meas(m_ori)// right left
1088 +//Symmetry
1089 alpha_r * dN_rl * ( var1(m_space.right() ,m_def.right() ) * var1(m_space.left(),m_def.left()).tr() ) * meas(m_ori)// right left
1090 ,
1091 alpha_r * dN_lr * var2( m_space.right(),m_space.right(),m_def.right(),usn(m_def.left() ).tr() ) * meas(m_ori) // right right
1092 +//Symmetry
1093 alpha_r * dN_rl * var2(m_space.right() ,m_space.right() ,m_def.right() ,usn(m_def.left()).tr() ) * meas(m_ori) // right right
1094 );
1095
1096 // Penalty of out-of-plane coupling
1097 // dW^pr / du_r --> second line
1098 m_assembler.assembleIfc(m_weakC1,
1099 alpha_r * ( var1(m_space.left(),m_def.left()) * usn(m_def.right()) ) * ( var1(m_space.left(),m_def.left()) * usn(m_def.right()) ).tr() * meas(m_ori) // left left
1100 ,
1101 alpha_r * ( var1(m_space.left(),m_def.left()) * usn(m_def.right()) ) * ( var1(m_space.right(),m_def.right()) * usn(m_def.left()) ).tr() * meas(m_ori) // left right
1102 ,
1103 alpha_r * ( var1(m_space.right(),m_def.right()) * usn(m_def.left()) ) * ( var1(m_space.left(),m_def.left()) * usn(m_def.right()) ).tr() * meas(m_ori) // right left
1104 ,
1105 alpha_r * ( var1(m_space.right(),m_def.right()) * usn(m_def.left()) ) * ( var1(m_space.right(),m_def.right()) * usn(m_def.left()) ).tr() * meas(m_ori) // right right
1106 ,
1107 // Symmetry
1108 alpha_r * ( var1(m_space.right(),m_def.right()) * usn(m_def.left()) ) * ( var1(m_space.right(),m_def.right()) * usn(m_def.left()) ).tr() * meas(m_ori) // right right
1109 ,
1110 alpha_r * ( var1(m_space.right(),m_def.right()) * usn(m_def.left()) ) * ( var1(m_space.left(),m_def.left()) * usn(m_def.right()) ).tr() * meas(m_ori) // right left
1111 ,
1112 alpha_r * ( var1(m_space.left(),m_def.left()) * usn(m_def.right()) ) * ( var1(m_space.right(),m_def.right()) * usn(m_def.left()) ).tr() * meas(m_ori) // left right
1113 ,
1114 alpha_r * ( var1(m_space.left(),m_def.left()) * usn(m_def.right()) ) * ( var1(m_space.left(),m_def.left()) * usn(m_def.right()) ).tr() * meas(m_ori) // left left
1115 );
1116
1117 // Penalty of in-plane coupling
1118 // dW^pr / du_r --> third line
1119 m_assembler.assembleIfc(m_weakC1,
1120 alpha_r * dnN_lr * ovar2(m_space.left(),m_space.left(),m_def.left(),usn(m_def.right()).tr()) * meas(m_ori) // left left
1121 + // Symmetry
1122 alpha_r * dnN_rl * ovar2(m_space.left(),m_space.left(),m_def.left(),usn(m_def.right()).tr()) * meas(m_ori) // left left
1123 ,
1124 alpha_r * dnN_lr * ( ovar1(m_space.left() ,m_def.left() ) * var1(m_space.right(),m_def.right()).tr() ) * meas(m_ori) // left right
1125 + // Symmetry
1126 alpha_r * dnN_rl * ( ovar1(m_space.left() ,m_def.left() ) * var1(m_space.right(),m_def.right()).tr() ) * meas(m_ori) // right left
1127 ,
1128 alpha_r * dnN_lr * ( ovar1(m_space.right(),m_def.right()) * var1(m_space.left() ,m_def.left() ).tr() ) * meas(m_ori) // right left
1129 + // Symmetry
1130 alpha_r * dnN_rl * ( ovar1(m_space.right(),m_def.right()) * var1(m_space.left() ,m_def.left() ).tr() ) * meas(m_ori) // right left
1131 ,
1132 alpha_r * dnN_lr * ovar2(m_space.right(),m_space.right(),m_def.right(),usn(m_def.left()).tr()) * meas(m_ori) // right right
1133 + // Symmetry
1134 alpha_r * dnN_rl * ovar2(m_space.right(),m_space.right(),m_def.right(),usn(m_def.left()).tr()) * meas(m_ori) // right right
1135 );
1136
1137 // Penalty of in-plane coupling
1138 // dW^pr / du_r --> fourth line
1139 m_assembler.assembleIfc(m_weakC1,
1140 alpha_r * ( ovar1(m_space.left() ,m_def.left() ) * usn(m_def.right()) ) * ( ovar1(m_space.left() ,m_def.left() ) * usn(m_def.right()) ).tr() * meas(m_ori) // left left
1141 + // Symmetry
1142 alpha_r * ( var1(m_space.left() ,m_def.left() ) * unv(m_def.right()) ) * ( var1(m_space.left() ,m_def.left() ) * unv(m_def.right()) ).tr() * meas(m_ori) // left left
1143 ,
1144 alpha_r * ( ovar1(m_space.left() ,m_def.left() ) * usn(m_def.right()) ) * ( var1(m_space.right(),m_def.right()) * unv(m_def.left() ) ).tr() * meas(m_ori) // left right
1145 + // Symmetry
1146 alpha_r * ( var1(m_space.left() ,m_def.left() ) * unv(m_def.right()) ) * ( ovar1(m_space.right(),m_def.right()) * usn(m_def.left() ) ).tr() * meas(m_ori) // left right
1147 ,
1148 alpha_r * ( var1(m_space.right(),m_def.right()) * unv(m_def.left() ) ) * ( ovar1(m_space.left() ,m_def.left() ) * usn(m_def.right()) ).tr() * meas(m_ori) // right left
1149 + // Symmetry
1150 alpha_r * ( ovar1(m_space.right(),m_def.right()) * usn(m_def.left() ) ) * ( var1(m_space.left() ,m_def.left() ) * unv(m_def.right()) ).tr() * meas(m_ori) // right left
1151 ,
1152 alpha_r * ( var1(m_space.right(),m_def.right()) * unv(m_def.left() ) ) * ( var1(m_space.right(),m_def.right()) * unv(m_def.left() ) ).tr() * meas(m_ori) // right right
1153 + // Symmetry
1154 alpha_r * ( ovar1(m_space.right(),m_def.right()) * usn(m_def.left() ) ) * ( ovar1(m_space.right(),m_def.right()) * usn(m_def.left() ) ).tr() * meas(m_ori) // right right
1155 );
1156}
1157
1158template <short_t d, class T, bool bending>
1159template <short_t _d, bool _matrix>
1160typename std::enable_if<(_d==3) && !_matrix, void>::type
1161gsThinShellAssembler<d, T, bending>::_assembleWeakIfc_impl(const gsFunctionSet<T> & deformed)
1162{
1163 geometryMap m_ori = m_assembler.getMap(m_patches);
1164 geometryMap m_def = m_assembler.getMap(deformed);
1165
1166 space m_space = m_assembler.trialSpace(0); // last argument is the space ID
1167
1168 gsMaterialMatrixIntegrate<T,MaterialOutput::MatrixA> m_mmA(m_materialMatrices,&m_patches,&deformed);
1169 gsMaterialMatrixIntegrate<T,MaterialOutput::MatrixD> m_mmD(m_materialMatrices,&m_patches,&deformed);
1170 auto mmA = m_assembler.getCoeff(m_mmA);
1171 auto mmD = m_assembler.getCoeff(m_mmD);
1172
1173 auto cart2cov = cartcov(m_ori);
1174 auto con2cartI = cartcon(m_ori).inv();
1175
1176 auto mmAcart = (con2cartI * reshape(mmA,3,3) * cart2cov);
1177 auto mmDcart = (con2cartI * reshape(mmD,3,3) * cart2cov);
1178
1179 element el = m_assembler.getElement();
1180 auto h = (el.area(m_ori.left()) + el.area(m_ori.right())) / 2;
1181 auto alpha_d = m_alpha_d_ifc * reshape(mmAcart,9,1).max().val() / h;
1182 auto alpha_r = m_alpha_r_ifc * reshape(mmDcart,9,1).max().val() / h;
1183
1184 auto du = ((m_def.left()-m_ori.left()) - (m_def.right()-m_ori.right()));
1185
1186 auto dN_lr = (usn(m_def.left()).tr()*usn(m_def.right())
1187 - usn(m_ori.left()).tr()*usn(m_ori.right())).val();
1188
1189 auto dN_rl = (usn(m_def.right()).tr()*usn(m_def.left())
1190 - usn(m_ori.right()).tr()*usn(m_ori.left())).val();
1191
1192 auto dnN_lr= (unv(m_def.left()).tr()*usn(m_def.right())
1193 - unv(m_ori.left()).tr()*usn(m_ori.right())).val();
1194
1195 auto dnN_rl= (unv(m_def.right()).tr()*usn(m_def.left())
1196 - unv(m_ori.right()).tr()*usn(m_ori.left())).val();
1197
1198 // C^0 coupling
1199 m_assembler.assembleIfc(m_weakC0,
1200 -alpha_d * m_space.left() * du * meas(m_ori)
1201 ,
1202 alpha_d * m_space.right()* du * meas(m_ori)
1203 );
1204
1205 // C^1 coupling
1206 m_assembler.assembleIfc(m_weakC1,
1207 -alpha_r * dN_lr * var1(m_space.left(),m_def.left()) * usn(m_def.right()) * meas(m_ori)
1208 ,
1209 -alpha_r * dN_lr * var1(m_space.right(),m_def.right()) * usn(m_def.left() ) * meas(m_ori)
1210 ,
1211 // Symmetry
1212 -alpha_r * dN_rl * var1(m_space.right(),m_def.right()) * usn(m_def.left()) * meas(m_ori)
1213 ,
1214 -alpha_r * dN_rl * var1(m_space.left(),m_def.left()) * usn(m_def.right() ) * meas(m_ori)
1215 );
1216
1217 // Penalty of in-plane coupling
1218 // dW^pr / du_r --> second line
1219 m_assembler.assembleIfc(m_weakC1,
1220 -alpha_r * dnN_lr* ovar1(m_space.left(),m_def.left()) * usn(m_def.right()) * meas(m_ori)
1221 ,
1222 -alpha_r * dnN_lr* var1(m_space.right(),m_def.right()) * unv(m_def.left() ) * meas(m_ori)
1223 ,
1224 // Symmetry
1225 -alpha_r * dnN_rl* ovar1(m_space.right(),m_def.right()) * usn(m_def.left()) * meas(m_ori)
1226 ,
1227 -alpha_r * dnN_rl* var1(m_space.left(),m_def.left()) * unv(m_def.right() ) * meas(m_ori)
1228 );
1229}
1230
1231template <short_t d, class T, bool bending>
1232template <short_t _d, bool _matrix>
1233typename std::enable_if<!(_d==3) && _matrix, void>::type
1234gsThinShellAssembler<d, T, bending>::_assembleWeakIfc_impl(const gsFunctionSet<T> & deformed)
1235{
1236 geometryMap m_ori = m_assembler.getMap(m_patches);
1237 geometryMap m_def = m_assembler.getMap(deformed);
1238
1239 space m_space = m_assembler.trialSpace(0); // last argument is the space ID
1240
1241 gsMaterialMatrixIntegrate<T,MaterialOutput::MatrixA> m_mmA(m_materialMatrices,&m_patches,&deformed);
1242 auto mmA = m_assembler.getCoeff(m_mmA);
1243
1244 auto cart2cov = cartcov(m_ori);
1245 auto con2cartI = cartcon(m_ori).inv();
1246
1247 auto mmAcart = (con2cartI * reshape(mmA,3,3) * cart2cov);
1248
1249 element el = m_assembler.getElement();
1250 auto h = (el.area(m_ori.left()) + el.area(m_ori.right())) / 2;
1251 auto alpha_d = m_alpha_d_ifc * reshape(mmAcart,9,1).max().val() / h;
1252
1253 auto du = ((m_def.left()-m_ori.left()) - (m_def.right()-m_ori.right()));
1254
1255 // C^0 coupling
1256 m_assembler.assembleIfc(m_weakC0,
1257 alpha_d * m_space.left() * m_space.left().tr() * meas(m_ori)
1258 ,
1259 -alpha_d * m_space.right()* m_space.left() .tr() * meas(m_ori)
1260 ,
1261 -alpha_d * m_space.left() * m_space.right().tr() * meas(m_ori)
1262 ,
1263 alpha_d * m_space.right()* m_space.right().tr() * meas(m_ori)
1264 );
1265
1266 // C^1 coupling DOES NOT CONTRIBUTE IN 2D PROBLEMS
1267}
1268
1269template <short_t d, class T, bool bending>
1270template <short_t _d, bool _matrix>
1271typename std::enable_if<!(_d==3) && !_matrix, void>::type
1272gsThinShellAssembler<d, T, bending>::_assembleWeakIfc_impl(const gsFunctionSet<T> & deformed)
1273{
1274 geometryMap m_ori = m_assembler.getMap(m_patches);
1275 geometryMap m_def = m_assembler.getMap(deformed);
1276
1277 space m_space = m_assembler.trialSpace(0); // last argument is the space ID
1278
1279 gsMaterialMatrixIntegrate<T,MaterialOutput::MatrixA> m_mmA(m_materialMatrices,&m_patches,&deformed);
1280 auto mmA = m_assembler.getCoeff(m_mmA);
1281
1282 auto cart2cov = cartcov(m_ori);
1283 auto con2cartI = cartcon(m_ori).inv();
1284
1285 auto mmAcart = (con2cartI * reshape(mmA,3,3) * cart2cov);
1286
1287 element el = m_assembler.getElement();
1288 auto h = (el.area(m_ori.left()) + el.area(m_ori.right())) / 2;
1289 auto alpha_d = m_alpha_d_ifc * reshape(mmAcart,9,1).max().val() / h;
1290
1291 auto du = ((m_def.left()-m_ori.left()) - (m_def.right()-m_ori.right()));
1292
1293 // C^0 coupling
1294 m_assembler.assembleIfc(m_weakC0,
1295 alpha_d * m_space.left() * du * meas(m_ori)
1296 ,
1297 -alpha_d * m_space.right()* du * meas(m_ori)
1298 );
1299
1300 // C^1 coupling DOES NOT CONTRIBUTE IN 2D PROBLEMS
1301}
1302
1303template<short_t d, class T, bool bending>
1304void gsThinShellAssembler<d, T, bending>::_assembleDirichlet()
1305{
1306 this->_getOptions();
1307 space m_space = m_assembler.trialSpace(0); // last argument is the space ID
1308 // if statement
1309 m_space.setup(m_bcs, dirichlet::l2Projection, m_continuity);
1310 // m_assembler.initSystem();
1311}
1312
1313template<short_t d, class T, bool bending>
1314void gsThinShellAssembler<d, T, bending>::homogenizeDirichlet()
1315{
1316 this->_getOptions();
1317 space m_space = m_assembler.trialSpace(0); // last argument is the space ID
1318 m_space.setup(m_bcs, dirichlet::homogeneous, m_continuity);
1319 // space m_space = m_assembler.trialSpace(0); // last argument is the space ID
1320 // const_cast<expr::gsFeSpace<T> & >(m_space).fixedPart().setZero();
1321}
1322
1323template<short_t d, class T, bool bending>
1324void gsThinShellAssembler<d, T, bending>::_applyLoads()
1325{
1326 gsMatrix<T> bVals;
1327 gsMatrix<index_t> acts,globalActs;
1328
1329 space m_space = m_assembler.trialSpace(0);
1330 m_mapper = m_space.mapper();
1331
1332 for (size_t i = 0; i< m_pLoads.numLoads(); ++i )
1333 {
1334 GISMO_ASSERT(m_pLoads[i].value.size()==d,"Point load has wrong dimension "<<m_pLoads[i].value.size()<<" instead of "<<d<<"\n");
1335 GISMO_ASSERT((size_t)m_pLoads[i].patch<m_patches.nPatches(),"Point load is defined on a patch with index "<<m_pLoads[i].patch<<" while the geometry has "<<m_patches.nPatches()<<" patches\n");
1336 // Compute actives and values of basis functions on point load location.
1337 if ( m_pLoads[i].parametric ) // in parametric space
1338 {
1339 if (const gsMappedBasis<2,T> * mbasis = dynamic_cast<const gsMappedBasis<2,T> * >(m_spaceBasis))
1340 {
1341 mbasis->active_into(m_pLoads[i].patch,m_pLoads[i].point, acts );
1342 mbasis->eval_into (m_pLoads[i].patch,m_pLoads[i].point, bVals );
1343 }
1344 else if (const gsMultiBasis<T> * mbasis = dynamic_cast<const gsMultiBasis<T> * >(m_spaceBasis))
1345 {
1346 mbasis->basis(m_pLoads[i].patch).active_into( m_pLoads[i].point, acts);
1347 mbasis->basis(m_pLoads[i].patch).eval_into ( m_pLoads[i].point, bVals);
1348 }
1349 else
1350 GISMO_ERROR("Basis type not understood");
1351 }
1352 else // in physical space
1353 {
1354 gsMatrix<T> forcePoint;
1355 m_patches.patch(m_pLoads[i].patch).invertPoints(m_pLoads[i].point,forcePoint);
1356
1357 if (const gsMappedBasis<2,T> * mbasis = dynamic_cast<const gsMappedBasis<2,T> * >(m_spaceBasis))
1358 {
1359 mbasis->active_into(m_pLoads[i].patch,forcePoint, acts );
1360 mbasis->eval_into (m_pLoads[i].patch,forcePoint, bVals );
1361 }
1362 else if (const gsMultiBasis<T> * mbasis = dynamic_cast<const gsMultiBasis<T> * >(m_spaceBasis))
1363 {
1364 mbasis->basis(m_pLoads[i].patch).active_into( forcePoint, acts);
1365 mbasis->basis(m_pLoads[i].patch).eval_into ( forcePoint, bVals);
1366 }
1367 else
1368 GISMO_ERROR("Basis type not understood");
1369 }
1370
1371 // Add the point load values in the right entries in the global RHS
1372 for (size_t j = 0; j< d; ++j)
1373 {
1374 if (m_pLoads[i].value[j] != 0.0)
1375 {
1376 m_mapper.localToGlobal(acts, m_pLoads[i].patch, globalActs,j);
1377 for (index_t k=0; k < globalActs.rows(); ++k)
1378 {
1379 if (m_mapper.is_free_index(globalActs(k,0)))
1380 m_rhs(globalActs(k,0), 0) += bVals(k,0) * m_pLoads[i].value[j];
1381 }
1382 }
1383 }
1384 }
1385}
1386
1387template<short_t d, class T, bool bending>
1388void gsThinShellAssembler<d, T, bending>::_applyMass()
1389{
1390 gsMatrix<T> bVals;
1391 gsMatrix<index_t> acts,globalActs;
1392
1393 space m_space = m_assembler.trialSpace(0);
1394 m_mapper = m_space.mapper();
1395
1396 GISMO_ASSERT(m_mass.rows()!=0,"Mass matrix must be assembled first");
1397
1398 for (size_t i = 0; i< m_pMass.numLoads(); ++i )
1399 {
1400 GISMO_ASSERT(m_pMass[i].value.size()==1,"Mass should be one-dimensional");
1401
1402 // Compute actives and values of basis functions on point load location.
1403 if ( m_pMass[i].parametric ) // in parametric space
1404 {
1405 m_basis.front().basis(m_pMass[i].patch).active_into( m_pMass[i].point, acts );
1406 m_basis.front().basis(m_pMass[i].patch).eval_into ( m_pMass[i].point, bVals);
1407 }
1408 else // in physical space
1409 {
1410 gsMatrix<T> forcePoint;
1411 m_patches.patch(m_pMass[i].patch).invertPoints(m_pMass[i].point,forcePoint);
1412 m_basis.front().basis(m_pMass[i].patch).active_into( forcePoint, acts );
1413 m_basis.front().basis(m_pMass[i].patch).eval_into ( forcePoint, bVals);
1414 }
1415
1416 // Add the point load values in the right entries in the global RHS
1417 for (size_t j = 0; j< d; ++j)
1418 {
1419 if (m_pMass[i].value[0] != 0.0)
1420 {
1421 m_mapper.localToGlobal(acts, m_pMass[i].patch, globalActs,j);
1422 for (index_t k=0; k < globalActs.rows(); ++k)
1423 {
1424 for (index_t l=0; l < globalActs.rows(); ++l)
1425 {
1426 if (m_mapper.is_free_index(globalActs(k,0)) && m_mapper.is_free_index(globalActs(l,0)))
1427 m_mass(globalActs(k,0), globalActs(l,0)) += bVals(k,0) * bVals(l,0) * m_pMass[i].value[0];
1428 }
1429 }
1430 }
1431 }
1432 }
1433}
1434
1435template<short_t d, class T, bool bending>
1436ThinShellAssemblerStatus gsThinShellAssembler<d, T, bending>::assembleMass(const bool lumped)
1437{
1438 m_assembler.cleanUp();
1439 GISMO_ENSURE(m_options.hasGroup("ExprAssembler"),"The option list does not contain options with the label 'ExprAssembler'!");
1440 m_assembler.setOptions(m_options.getGroup("ExprAssembler"));
1441
1442 geometryMap m_ori = m_assembler.getMap(m_patches);
1443
1444 // Initialize stystem
1445 m_assembler.initSystem();
1446
1447 gsMaterialMatrixIntegrate<T,MaterialOutput::Density> m_mm(m_materialMatrices,&m_patches);
1448 auto mm0 = m_assembler.getCoeff(m_mm);
1449
1450 space m_space = m_assembler.trialSpace(0);
1451 m_space.setup(m_bcs, dirichlet::homogeneous, m_continuity);
1452 // this->homogenizeDirichlet();
1453
1454 gsExprEvaluator<T> ev(m_assembler);
1455 gsVector<> pt(2);
1456 pt.setConstant(.25);
1457
1458 try
1459 {
1460 // assemble system
1461 if (!lumped)
1462 m_assembler.assemble(mm0.val()*m_space*m_space.tr()*meas(m_ori));
1463 else
1464 m_assembler.assemble(mm0.val()*(m_space.rowSum())*meas(m_ori));
1465
1466 m_mass = m_assembler.matrix();
1467
1468/* // assemble system
1469 if (!lumped)
1470 {
1471 m_assembler.assemble(mm0.val()*m_space*m_space.tr()*meas(m_ori));
1472 m_mass = m_assembler.matrix();
1473 this->_applyMass();
1474 }
1475 else
1476 {
1477 // To do: add point masses in lumped case
1478 m_assembler.assemble(mm0.val()*(m_space.rowSum())*meas(m_ori));
1479 m_rhs = m_assembler.rhs();
1480 }
1481 m_mass = m_assembler.matrix();
1482
1483 this->_applyMass();
1484 m_status = ThinShellAssemblerStatus::Success;*/
1485 }
1486 catch (...)
1487 {
1488 m_assembler.cleanUp();
1489 m_status = ThinShellAssemblerStatus::AssemblyError;
1490 }
1491 return m_status;
1492}
1493
1494// legacy
1495template<short_t d, class T, bool bending>
1496ThinShellAssemblerStatus gsThinShellAssembler<d, T, bending>::assembleFoundation()
1497{
1498 m_assembler.cleanUp();
1499 GISMO_ENSURE(m_options.hasGroup("ExprAssembler"),"The option list does not contain options with the label 'ExprAssembler'!");
1500 m_assembler.setOptions(m_options.getGroup("ExprAssembler"));
1501
1502 geometryMap m_ori = m_assembler.getMap(m_patches);
1503
1504 // Initialize stystem
1505 m_assembler.initSystem();
1506 auto m_foundation = m_assembler.getCoeff(*m_foundFun, m_ori);
1507 GISMO_ASSERT(m_foundFun->targetDim()==3,"Foundation function has dimension "<<m_foundFun->targetDim()<<", but expected 3");
1508
1509 space m_space = m_assembler.trialSpace(0);
1510
1511 try
1512 {
1513 m_assembler.assemble(m_space * m_foundation.asDiag() * m_space.tr() * meas(m_ori));
1514 m_status = ThinShellAssemblerStatus::Success;
1515 }
1516 catch (...)
1517 {
1518 m_assembler.cleanUp();
1519 m_status = ThinShellAssemblerStatus::AssemblyError;
1520 }
1521 return m_status;
1522}
1523
1524template<short_t d, class T, bool bending>
1525ThinShellAssemblerStatus gsThinShellAssembler<d, T, bending>::assemble()
1526{
1527 return assemble_impl<d, bending>();
1528}
1529
1537template <short_t d, typename T, bool bending>
1538template <short_t _d, bool _bending>
1539typename std::enable_if<(_d==3) && _bending, ThinShellAssemblerStatus>::type
1540gsThinShellAssembler<d, T, bending>::assemble_impl()
1541{
1542 this->_getOptions();
1543
1544 m_assembler.cleanUp();
1545 GISMO_ENSURE(m_options.hasGroup("ExprAssembler"),"The option list does not contain options with the label 'ExprAssembler'!");
1546 m_assembler.setOptions(m_options.getGroup("ExprAssembler"));
1547
1548 // Linear assembly: deformed and undeformed geometries are the same
1549 gsMultiPatch<T> & defpatches = m_patches;
1550 geometryMap m_ori = m_assembler.getMap(m_patches);
1551 geometryMap m_def = m_assembler.getMap(defpatches);
1552
1553 // Initialize stystem
1554 m_assembler.initSystem();
1555 m_assembler.initVector(1);
1556
1557 gsMaterialMatrixIntegrate<T,MaterialOutput::MatrixA> m_mmA(m_materialMatrices,&m_patches,&defpatches);
1558 gsMaterialMatrixIntegrate<T,MaterialOutput::MatrixB> m_mmB(m_materialMatrices,&m_patches,&defpatches);
1559 gsMaterialMatrixIntegrate<T,MaterialOutput::MatrixC> m_mmC(m_materialMatrices,&m_patches,&defpatches);
1560 gsMaterialMatrixIntegrate<T,MaterialOutput::MatrixD> m_mmD(m_materialMatrices,&m_patches,&defpatches);
1561 auto mmA = m_assembler.getCoeff(m_mmA);
1562 auto mmB = m_assembler.getCoeff(m_mmB);
1563 auto mmC = m_assembler.getCoeff(m_mmC);
1564 auto mmD = m_assembler.getCoeff(m_mmD);
1565
1566 gsFunctionExpr<> mult2t("1","0","0","0","1","0","0","0","2",2);
1567 auto m_m2 = m_assembler.getCoeff(mult2t);
1568
1569 space m_space = m_assembler.trialSpace(0);
1570
1571 auto m_physforce = m_assembler.getCoeff(*m_forceFun,m_ori); // force defined in physical domain
1572 auto m_parforce = m_assembler.getCoeff(*m_forceFun); // force defined in parametric domain
1573
1574 auto m_Em_der = flat( jac(m_def).tr() * jac(m_space) );
1575 auto m_Ef_der = -( deriv2(m_space,sn(m_def).normalized().tr() ) + deriv2(m_def,var1(m_space,m_def) ) ) * reshape(m_m2,3,3);
1576
1577 auto m_N_der = m_Em_der * reshape(mmA,3,3) + m_Ef_der * reshape(mmB,3,3);
1578 auto m_M_der = m_Em_der * reshape(mmC,3,3) + m_Ef_der * reshape(mmD,3,3);
1579
1580 try
1581 {
1582 if (m_foundInd)
1583 {
1584 this->_assembleFoundation<true>(*m_foundFun);
1585 this->_assembleFoundation<false>(*m_foundFun);
1586 }
1587 if (m_pressInd)
1588 {
1589 this->_assemblePressure<true>(*m_pressFun);
1590 this->_assemblePressure<false>(*m_pressFun);
1591 }
1592
1593 m_assembler.assemble(
1594 (
1595 m_N_der * m_Em_der.tr()
1596 +
1597 m_M_der * m_Ef_der.tr()
1598 ) * meas(m_ori)
1599 );
1600
1601 if (m_parametricForce) // Assemble the force defined in the parameter domain
1602 m_assembler.assemble(m_space * m_parforce * meas(m_ori));
1603 else // Assemble the force defined in the physical domain
1604 m_assembler.assemble(m_space * m_physforce * meas(m_ori));
1605
1606 this->_assembleWeakBCs<true>();
1607 this->_assembleWeakBCs<false>();
1608 this->_assembleWeakIfc<true>();
1609 this->_assembleWeakIfc<false>();
1610 this->_assembleNeumann();
1611
1612 // Assemble the loads
1613 if ( m_pLoads.numLoads() != 0 )
1614 {
1615 m_rhs = m_assembler.rhs();
1616 _applyLoads();
1617 }
1618 m_status = ThinShellAssemblerStatus::Success;
1619 }
1620 catch (...)
1621 {
1622 m_assembler.cleanUp();
1623 m_status = ThinShellAssemblerStatus::AssemblyError;
1624 }
1625 return m_status;
1626}
1627
1628template <short_t d, typename T, bool bending>
1629template <short_t _d, bool _bending>
1630typename std::enable_if<!(_d==3 && _bending), ThinShellAssemblerStatus>::type
1631gsThinShellAssembler<d, T, bending>::assemble_impl()
1632{
1633 this->_getOptions();
1634
1635 m_assembler.cleanUp();
1636 GISMO_ENSURE(m_options.hasGroup("ExprAssembler"),"The option list does not contain options with the label 'ExprAssembler'!");
1637 m_assembler.setOptions(m_options.getGroup("ExprAssembler"));
1638
1639 // Linear assembly: deformed and undeformed geometries are the same
1640 gsMultiPatch<T> & defpatches = m_patches;
1641 geometryMap m_ori = m_assembler.getMap(m_patches);
1642 geometryMap m_def = m_assembler.getMap(defpatches);
1643
1644 // Initialize stystem
1645 m_assembler.initSystem();
1646
1647 gsMaterialMatrixIntegrate<T,MaterialOutput::MatrixA> m_mmA(m_materialMatrices,&m_patches,&defpatches);
1648 auto mmA = m_assembler.getCoeff(m_mmA);
1649
1650 space m_space = m_assembler.trialSpace(0);
1651 auto m_physforce = m_assembler.getCoeff(*m_forceFun,m_ori); // force defined in physical domain
1652 auto m_parforce = m_assembler.getCoeff(*m_forceFun); // force defined in parametric domain
1653
1654 auto jacG = jac(m_def);
1655 auto m_Em_der = flat( jacG.tr() * jac(m_space) ) ; //[checked]
1656 auto m_N_der = m_Em_der * reshape(mmA,3,3);
1657
1658 try
1659 {
1660 if (m_foundInd)
1661 {
1662 this->_assembleFoundation<true>(*m_foundFun);
1663 this->_assembleFoundation<false>(*m_foundFun);
1664 }
1665 if (m_pressInd)
1666 {
1667 this->_assemblePressure<true>(*m_pressFun);
1668 this->_assemblePressure<false>(*m_pressFun);
1669 }
1670
1671 m_assembler.assemble(
1672 (
1673 m_N_der * m_Em_der.tr()
1674 ) * meas(m_ori)
1675 );
1676
1677 if (m_parametricForce) // Assemble the force defined in the parameter domain
1678 m_assembler.assemble(m_space * m_parforce * meas(m_ori));
1679 else // Assemble the force defined in the physical domain
1680 m_assembler.assemble(m_space * m_physforce * meas(m_ori));
1681
1682 this->_assembleWeakBCs<true>();
1683 this->_assembleWeakBCs<false>();
1684 this->_assembleWeakIfc<true>();
1685 this->_assembleWeakIfc<false>();
1686 this->_assembleNeumann();
1687
1688 // Assemble the loads
1689 if ( m_pLoads.numLoads() != 0 )
1690 {
1691 m_rhs = m_assembler.rhs();
1692 _applyLoads();
1693 }
1694
1695 m_status = ThinShellAssemblerStatus::Success;
1696 }
1697 catch (...)
1698 {
1699 m_assembler.cleanUp();
1700 m_status = ThinShellAssemblerStatus::AssemblyError;
1701 }
1702 return m_status;
1703}
1704
1705template<short_t d, class T, bool bending>
1706ThinShellAssemblerStatus gsThinShellAssembler<d, T, bending>::assembleMatrix(const gsFunctionSet<T> & deformed)
1707{
1708 return assembleMatrix_impl<d, bending>(deformed);
1709}
1710
1711template <short_t d, typename T, bool bending>
1712template <short_t _d, bool _bending>
1713typename std::enable_if<(_d==3) && _bending, ThinShellAssemblerStatus>::type
1714gsThinShellAssembler<d, T, bending>::assembleMatrix_impl(const gsFunctionSet<T> & deformed)
1715{
1716 m_assembler.cleanUp();
1717 GISMO_ENSURE(m_options.hasGroup("ExprAssembler"),"The option list does not contain options with the label 'ExprAssembler'!");
1718 m_assembler.setOptions(m_options.getGroup("ExprAssembler"));
1719
1720 geometryMap m_ori = m_assembler.getMap(m_patches);
1721 geometryMap m_def = m_assembler.getMap(deformed);
1722
1723 // Initialize matrix
1724 m_assembler.initSystem();
1725 m_assembler.initMatrix();
1726
1727 gsMaterialMatrixIntegrate<T,MaterialOutput::MatrixA> m_mmA(m_materialMatrices,&m_patches,&deformed);
1728 gsMaterialMatrixIntegrate<T,MaterialOutput::MatrixB> m_mmB(m_materialMatrices,&m_patches,&deformed);
1729 gsMaterialMatrixIntegrate<T,MaterialOutput::MatrixC> m_mmC(m_materialMatrices,&m_patches,&deformed);
1730 gsMaterialMatrixIntegrate<T,MaterialOutput::MatrixD> m_mmD(m_materialMatrices,&m_patches,&deformed);
1731 gsMaterialMatrixIntegrate<T,MaterialOutput::VectorN> m_S0(m_materialMatrices,&m_patches,&deformed);
1732 gsMaterialMatrixIntegrate<T,MaterialOutput::VectorM> m_S1(m_materialMatrices,&m_patches,&deformed);
1733 auto mmA = m_assembler.getCoeff(m_mmA);
1734 auto mmB = m_assembler.getCoeff(m_mmB);
1735 auto mmC = m_assembler.getCoeff(m_mmC);
1736 auto mmD = m_assembler.getCoeff(m_mmD);
1737 auto S0 = m_assembler.getCoeff(m_S0);
1738 auto S1 = m_assembler.getCoeff(m_S1);
1739
1740 gsFunctionExpr<T> mult2t("1","0","0","0","1","0","0","0","2",2);
1741 auto m_m2 = m_assembler.getCoeff(mult2t);
1742
1743 space m_space = m_assembler.trialSpace(0);
1744
1745 this->homogenizeDirichlet();
1746
1747 auto m_N = S0.tr();
1748 auto m_Em_der = flat( jac(m_def).tr() * jac(m_space) ) ; //[checked]
1749 auto m_Em_der2 = flatdot( jac(m_space),jac(m_space).tr(), m_N ); //[checked]
1750
1751
1752 auto m_M = S1.tr(); // output is a column
1753 auto m_Ef_der = -( deriv2(m_space,sn(m_def).normalized().tr() ) + deriv2(m_def,var1(m_space,m_def) ) ) * reshape(m_m2,3,3); //[checked]
1754 auto m_Ef_der2 = -(flatdot2( deriv2(m_space), var1(m_space,m_def).tr(), m_M ).symmetrize()
1755 + var2deriv2(m_space,m_space,m_def, m_M ));
1756
1757 auto m_N_der = m_Em_der * reshape(mmA,3,3) + m_Ef_der * reshape(mmB,3,3);
1758 auto m_M_der = m_Em_der * reshape(mmC,3,3) + m_Ef_der * reshape(mmD,3,3);
1759
1760 try
1761 {
1762 if (m_foundInd) this->_assembleFoundation<true>(*m_foundFun,deformed);
1763 if (m_pressInd) this->_assemblePressure<true>(*m_pressFun,deformed);
1764
1765 // Assemble matrix
1766 m_assembler.assemble(
1767 (
1768 m_N_der * m_Em_der.tr()
1769 +
1770 m_Em_der2
1771 +
1772 m_M_der * m_Ef_der.tr()
1773 +
1774 m_Ef_der2
1775 ) * meas(m_ori)
1776 );
1777
1778 this->_assembleWeakBCs<true>(deformed);
1779 this->_assembleWeakIfc<true>(deformed);
1780
1781 m_status = ThinShellAssemblerStatus::Success;
1782 }
1783 catch (...)
1784 {
1785 m_assembler.cleanUp();
1786 m_status = ThinShellAssemblerStatus::AssemblyError;
1787 }
1788 return m_status;
1789}
1790
1791template <short_t d, typename T, bool bending>
1792template <short_t _d, bool _bending>
1793typename std::enable_if<!(_d==3 && _bending), ThinShellAssemblerStatus>::type
1794gsThinShellAssembler<d, T, bending>::assembleMatrix_impl(const gsFunctionSet<T> & deformed)
1795{
1796 m_assembler.cleanUp();
1797 GISMO_ENSURE(m_options.hasGroup("ExprAssembler"),"The option list does not contain options with the label 'ExprAssembler'!");
1798 m_assembler.setOptions(m_options.getGroup("ExprAssembler"));
1799
1800 geometryMap m_ori = m_assembler.getMap(m_patches);
1801 geometryMap m_def = m_assembler.getMap(deformed);
1802
1803 // Initialize matrix
1804 m_assembler.initSystem();
1805 m_assembler.initMatrix();
1806
1807 gsMaterialMatrixIntegrate<T,MaterialOutput::MatrixA> m_mmA(m_materialMatrices,&m_patches,&deformed);
1808 gsMaterialMatrixIntegrate<T,MaterialOutput::VectorN> m_S0(m_materialMatrices,&m_patches,&deformed);
1809 auto mmA = m_assembler.getCoeff(m_mmA);
1810 auto S0 = m_assembler.getCoeff(m_S0);
1811
1812 space m_space = m_assembler.trialSpace(0);
1813
1814
1815 this->homogenizeDirichlet();
1816
1817 auto m_N = S0.tr();
1818 auto m_Em_der = flat( jac(m_def).tr() * jac(m_space) ) ; //[checked]
1819 auto m_Em_der2 = flatdot( jac(m_space),jac(m_space).tr(), m_N ); //[checked]
1820
1821 auto m_N_der = m_Em_der * reshape(mmA,3,3);
1822
1823 // Assemble matrix
1824 try
1825 {
1826 if (m_foundInd) this->_assembleFoundation<true>(*m_foundFun,deformed);
1827 if (m_pressInd) this->_assemblePressure<true>(*m_pressFun,deformed);
1828
1829 m_assembler.assemble(
1830 (
1831 m_N_der * m_Em_der.tr()
1832 +
1833 m_Em_der2
1834 ) * meas(m_ori)
1835 );
1836 this->_assembleWeakBCs<true>(deformed);
1837 this->_assembleWeakIfc<true>(deformed);
1838
1839 m_status = ThinShellAssemblerStatus::Success;
1840 }
1841 catch (...) // add specific cases?
1842 {
1843 m_assembler.cleanUp();
1844 m_status = ThinShellAssemblerStatus::AssemblyError;
1845 }
1846 return m_status;
1847}
1848
1849template<short_t d, class T, bool bending>
1850ThinShellAssemblerStatus gsThinShellAssembler<d, T, bending>::assembleMatrix(const gsMatrix<T> & solVector)
1851{
1852 gsMultiPatch<T> def;
1853 constructSolution(solVector, def);
1854 return assembleMatrix(def);
1855}
1856
1857template<short_t d, class T, bool bending>
1858ThinShellAssemblerStatus gsThinShellAssembler<d, T, bending>::assembleMatrix(const gsFunctionSet<T> & deformed, const gsFunctionSet<T> & previous, gsMatrix<T> & update)
1859{
1860 return assembleMatrix_impl<d, bending>(deformed, previous, update);
1861}
1862
1863template <short_t d, typename T, bool bending>
1864template <short_t _d, bool _bending>
1865typename std::enable_if<(_d==3) && _bending, ThinShellAssemblerStatus>::type
1866gsThinShellAssembler<d, T, bending>::assembleMatrix_impl(const gsFunctionSet<T> & deformed, const gsFunctionSet<T> & previous, gsMatrix<T> & update)
1867{
1868 m_assembler.cleanUp();
1869 GISMO_ENSURE(m_options.hasGroup("ExprAssembler"),"The option list does not contain options with the label 'ExprAssembler'!");
1870 m_assembler.setOptions(m_options.getGroup("ExprAssembler"));
1871
1872 geometryMap m_ori = m_assembler.getMap(m_patches);
1873 geometryMap m_def = m_assembler.getMap(deformed);
1874 geometryMap m_prev = m_assembler.getMap(previous);
1875 // Initialize matrix
1876 m_assembler.initSystem();
1877 m_assembler.initMatrix();
1878
1879 gsMaterialMatrixIntegrate<T,MaterialOutput::MatrixA> m_mmA(m_materialMatrices,&m_patches,&deformed);
1880 gsMaterialMatrixIntegrate<T,MaterialOutput::MatrixB> m_mmB(m_materialMatrices,&m_patches,&deformed);
1881 gsMaterialMatrixIntegrate<T,MaterialOutput::MatrixC> m_mmC(m_materialMatrices,&m_patches,&deformed);
1882 gsMaterialMatrixIntegrate<T,MaterialOutput::MatrixD> m_mmD(m_materialMatrices,&m_patches,&deformed);
1883 gsMaterialMatrixIntegrate<T,MaterialOutput::VectorN> m_S0(m_materialMatrices,&m_patches,&deformed);
1884 gsMaterialMatrixIntegrate<T,MaterialOutput::VectorM> m_S1(m_materialMatrices,&m_patches,&deformed);
1885 auto mmA = m_assembler.getCoeff(m_mmA);
1886 auto mmB = m_assembler.getCoeff(m_mmB);
1887 auto mmC = m_assembler.getCoeff(m_mmC);
1888 auto mmD = m_assembler.getCoeff(m_mmD);
1889 // auto S0 = m_assembler.getCoeff(m_S0);
1890 // auto S1 = m_assembler.getCoeff(m_S1);
1891
1892 gsMaterialMatrixIntegrate<T,MaterialOutput::MatrixA> m_mmAd(m_materialMatrices,&m_patches,&previous);
1893 gsMaterialMatrixIntegrate<T,MaterialOutput::MatrixB> m_mmBd(m_materialMatrices,&m_patches,&previous);
1894 gsMaterialMatrixIntegrate<T,MaterialOutput::MatrixC> m_mmCd(m_materialMatrices,&m_patches,&previous);
1895 gsMaterialMatrixIntegrate<T,MaterialOutput::MatrixD> m_mmDd(m_materialMatrices,&m_patches,&previous);
1896 gsMaterialMatrixIntegrate<T,MaterialOutput::VectorN> m_S0d(m_materialMatrices,&m_patches,&previous);
1897 gsMaterialMatrixIntegrate<T,MaterialOutput::VectorM> m_S1d(m_materialMatrices,&m_patches,&previous);
1898 auto mmAp = m_assembler.getCoeff(m_mmAd);
1899 auto mmBp = m_assembler.getCoeff(m_mmBd);
1900 auto mmCp = m_assembler.getCoeff(m_mmCd);
1901 auto mmDp = m_assembler.getCoeff(m_mmDd);
1902 auto S0 = m_assembler.getCoeff(m_S0d);
1903 auto S1 = m_assembler.getCoeff(m_S1d);
1904
1905 gsFunctionExpr<> mult2t("1","0","0","0","1","0","0","0","2",2);
1906 auto m_m2 = m_assembler.getCoeff(mult2t);
1907
1908 space m_space = m_assembler.trialSpace(0);
1909 solution m_du = m_assembler.getSolution(m_space,update);
1910
1911 this->homogenizeDirichlet();
1912
1913 auto m_E_mc = flat( jac(m_prev).tr() * grad(m_du) ) ; //[checked]
1914 auto m_E_fc = -( deriv2(m_du,sn(m_prev).normalized().tr() ) + deriv2(m_prev,var1(m_du,m_prev) ) ) * reshape(m_m2,3,3); //[checked]
1915 auto m_N_c = m_E_mc * reshape(mmAp,3,3) + m_E_fc * reshape(mmBp,3,3);
1916 auto m_M_c = m_E_mc * reshape(mmCp,3,3) + m_E_fc * reshape(mmDp,3,3);
1917
1918 auto m_N = S0.tr() + m_N_c;
1919 auto m_Em_der = flat( jac(m_def).tr() * jac(m_space) ) ; //[checked]
1920 auto m_Em_der2 = flatdot( jac(m_space),jac(m_space).tr(), m_N ); //[checked]
1921
1922 auto m_M = S1.tr() + m_M_c; // output is a column
1923 auto m_Ef_der = -( deriv2(m_space,sn(m_def).normalized().tr() ) + deriv2(m_def,var1(m_space,m_def) ) ) * reshape(m_m2,3,3); //[checked]
1924 auto m_Ef_der2 = -(flatdot2( deriv2(m_space), var1(m_space,m_def).tr(), m_M ).symmetrize()
1925 + var2deriv2(m_space,m_space,m_def, m_M ));
1926
1927 auto m_N_der = m_Em_der * reshape(mmA,3,3) + m_Ef_der * reshape(mmB,3,3);
1928 auto m_M_der = m_Em_der * reshape(mmC,3,3) + m_Ef_der * reshape(mmD,3,3);
1929
1930 try
1931 {
1932 if (m_foundInd) this->_assembleFoundation<true>(*m_foundFun,deformed);
1933 if (m_pressInd) this->_assemblePressure<true>(*m_pressFun,deformed);
1934
1935 // Assemble matrix
1936 m_assembler.assemble(
1937 (
1938 m_N_der * m_Em_der.tr()
1939 +
1940 m_Em_der2
1941 +
1942 m_M_der * m_Ef_der.tr()
1943 +
1944 m_Ef_der2
1945 ) * meas(m_ori)
1946 );
1947 this->_assembleWeakBCs<true>(deformed);
1948 this->_assembleWeakIfc<true>(deformed);
1949
1950 m_status = ThinShellAssemblerStatus::Success;
1951 }
1952 catch (...)
1953 {
1954 m_assembler.cleanUp();
1955 m_status = ThinShellAssemblerStatus::AssemblyError;
1956 }
1957 return m_status;
1958}
1959
1960// template <short_t d, typename T, bool bending>
1961// template <short_t _d, bool _bending>
1962// typename std::enable_if<!(_d==3 && _bending), ThinShellAssemblerStatus>::type
1963// gsThinShellAssembler<d, T, bending>::assembleMatrix_impl(const gsMultiPatch<T> & deformed, const gsMultiPatch<T> & previous, gsMatrix<T> & update)
1964// {
1965// GISMO_NO_IMPLEMENTATION;
1966// }
1967
1968template<short_t d, class T, bool bending>
1969ThinShellAssemblerStatus gsThinShellAssembler<d, T, bending>::assembleMatrix(const gsMatrix<T> & solVector, const gsMatrix<T> & prevVector)
1970{
1971 // gsMultiPatch<T> deformed;
1972 // constructSolution(solVector, deformed);
1973 // assembleMatrix(deformed);
1974
1975 gsMultiPatch<T> def, it;
1976 constructSolution(solVector, def);
1977 constructSolution(prevVector, it);
1978 gsMatrix<T> update = solVector - prevVector;
1979 return assembleMatrix(def,it,update);
1980}
1981
1982template<short_t d, class T, bool bending>
1983ThinShellAssemblerStatus gsThinShellAssembler<d, T, bending>::assembleVector(const gsFunctionSet<T> & deformed, const bool homogenize)
1984{
1985 return assembleVector_impl<d, bending>(deformed,homogenize);
1986}
1987
1988template <short_t d, typename T, bool bending>
1989template <short_t _d, bool _bending>
1990typename std::enable_if<(_d==3) && _bending, ThinShellAssemblerStatus>::type
1991gsThinShellAssembler<d, T, bending>::assembleVector_impl(const gsFunctionSet<T> & deformed, const bool homogenize)
1992{
1993 m_assembler.cleanUp();
1994 GISMO_ENSURE(m_options.hasGroup("ExprAssembler"),"The option list does not contain options with the label 'ExprAssembler'!");
1995 m_assembler.setOptions(m_options.getGroup("ExprAssembler"));
1996
1997 geometryMap m_ori = m_assembler.getMap(m_patches);
1998 geometryMap m_def = m_assembler.getMap(deformed);
1999
2000 // Initialize vector
2001 m_assembler.initVector(1);
2002
2003 gsMaterialMatrixIntegrate<T,MaterialOutput::VectorN> m_S0(m_materialMatrices,&m_patches,&deformed);
2004 gsMaterialMatrixIntegrate<T,MaterialOutput::VectorM> m_S1(m_materialMatrices,&m_patches,&deformed);
2005 auto S0 = m_assembler.getCoeff(m_S0);
2006 auto S1 = m_assembler.getCoeff(m_S1);
2007
2008 gsFunctionExpr<> mult2t("1","0","0","0","1","0","0","0","2",2);
2009 auto m_m2 = m_assembler.getCoeff(mult2t);
2010
2011 space m_space = m_assembler.trialSpace(0);
2012 auto m_physforce = m_assembler.getCoeff(*m_forceFun,m_ori); // force defined in physical domain
2013 auto m_parforce = m_assembler.getCoeff(*m_forceFun); // force defined in parametric domain
2014
2015
2016 if (homogenize) this->homogenizeDirichlet();
2017 else this->_assembleDirichlet();
2018
2019 auto m_N = S0.tr();
2020 auto m_Em_der = flat( jac(m_def).tr() * jac(m_space) ) ;
2021
2022 auto m_M = S1.tr(); // output is a column
2023 auto m_Ef_der = -( deriv2(m_space,sn(m_def).normalized().tr() ) + deriv2(m_def,var1(m_space,m_def) ) ) * reshape(m_m2,3,3); //[checked]
2024
2025 try
2026 {
2027 if (m_foundInd) this->_assembleFoundation<false>(*m_foundFun,deformed);
2028 if (m_pressInd) this->_assemblePressure<false>(*m_pressFun,deformed);
2029
2030 // Assemble vector
2032 if (m_parametricForce) // Assemble the force defined in the parameter domain
2033 m_assembler.assemble(m_space * m_parforce * meas(m_ori));
2034 else // Assemble the force defined in the physical domain
2035 m_assembler.assemble(m_space * m_physforce * meas(m_ori));
2036
2038 m_assembler.assemble(
2039 (
2040 - ( ( m_N * m_Em_der.tr() + m_M * m_Ef_der.tr() ) * meas(m_ori) ).tr()
2041 )
2042 );
2043
2044 this->_assembleWeakBCs<false>(deformed);
2045 this->_assembleWeakIfc<false>(deformed);
2046 this->_assembleNeumann();
2047
2048 // Assemble the loads
2049 if ( m_pLoads.numLoads() != 0 )
2050 {
2051 m_rhs = m_assembler.rhs();
2052 _applyLoads();
2053 }
2054
2055 m_status = ThinShellAssemblerStatus::Success;
2056 }
2057 catch (...)
2058 {
2059 m_assembler.cleanUp();
2060 m_status = ThinShellAssemblerStatus::AssemblyError;
2061 }
2062 return m_status;
2063}
2064
2065template <short_t d, typename T, bool bending>
2066template <short_t _d, bool _bending>
2067typename std::enable_if<!(_d==3 && _bending), ThinShellAssemblerStatus>::type
2068gsThinShellAssembler<d, T, bending>::assembleVector_impl(const gsFunctionSet<T> & deformed, const bool homogenize)
2069{
2070 m_assembler.cleanUp();
2071 GISMO_ENSURE(m_options.hasGroup("ExprAssembler"),"The option list does not contain options with the label 'ExprAssembler'!");
2072 m_assembler.setOptions(m_options.getGroup("ExprAssembler"));
2073
2074 geometryMap m_ori = m_assembler.getMap(m_patches);
2075 geometryMap m_def = m_assembler.getMap(deformed);
2076
2077 // Initialize vector
2078 m_assembler.initVector(1);
2079
2080 gsMaterialMatrixIntegrate<T,MaterialOutput::VectorN> m_S0(m_materialMatrices,&m_patches,&deformed);
2081 auto S0 = m_assembler.getCoeff(m_S0);
2082
2083 space m_space = m_assembler.trialSpace(0);
2084 auto m_physforce = m_assembler.getCoeff(*m_forceFun,m_ori); // force defined in physical domain
2085 auto m_parforce = m_assembler.getCoeff(*m_forceFun); // force defined in parametric domain
2086
2087 if (homogenize) this->homogenizeDirichlet();
2088 else this->_assembleDirichlet();
2089
2090 auto m_N = S0.tr();
2091 auto m_Em_der = flat( jac(m_def).tr() * jac(m_space) ) ;
2092
2093 try
2094 {
2095 if (m_foundInd) this->_assembleFoundation<false>(*m_foundFun,deformed);
2096 if (m_pressInd) this->_assemblePressure<false>(*m_pressFun,deformed);
2097
2098 // Assemble vector
2100 if (m_parametricForce) // Assemble the force defined in the parameter domain
2101 m_assembler.assemble(m_space * m_parforce * meas(m_ori));
2102 else // Assemble the force defined in the physical domain
2103 m_assembler.assemble(m_space * m_physforce * meas(m_ori));
2104
2106 m_assembler.assemble(
2107 (
2108 - ( ( m_N * m_Em_der.tr() ) * meas(m_ori) ).tr()
2109 )
2110 );
2111
2112 this->_assembleWeakBCs<false>(deformed);
2113 this->_assembleWeakIfc<false>(deformed);
2114 this->_assembleNeumann();
2115
2116 // Assemble the loads
2117 if ( m_pLoads.numLoads() != 0 )
2118 {
2119 m_rhs = m_assembler.rhs();
2120 _applyLoads();
2121 }
2122
2123 m_status = ThinShellAssemblerStatus::Success;
2124 }
2125 catch (...)
2126 {
2127 m_assembler.cleanUp();
2128 m_status = ThinShellAssemblerStatus::AssemblyError;
2129 }
2130 return m_status;
2131}
2132
2133template<short_t d, class T, bool bending>
2134ThinShellAssemblerStatus gsThinShellAssembler<d, T, bending>::assemblePressureMatrix(const gsFunction<T> & pressFun)
2135{
2136 try
2137 {
2138 this->_assemblePressure<true>(pressFun);
2139 m_status = ThinShellAssemblerStatus::Success;
2140 }
2141 catch (...)
2142 {
2143 m_status = ThinShellAssemblerStatus::AssemblyError;
2144 }
2145 return m_status;
2146}
2147
2148template<short_t d, class T, bool bending>
2149ThinShellAssemblerStatus gsThinShellAssembler<d, T, bending>::assemblePressureMatrix(const T pressure)
2150{
2151 gsConstantFunction<T> pressFun(pressure,d);
2152 try
2153 {
2154 this->_assemblePressure<true>(pressFun);
2155 m_status = ThinShellAssemblerStatus::Success;
2156 }
2157 catch (...)
2158 {
2159 m_status = ThinShellAssemblerStatus::AssemblyError;
2160 }
2161 return m_status;
2162}
2163
2164template<short_t d, class T, bool bending>
2165ThinShellAssemblerStatus gsThinShellAssembler<d, T, bending>::assemblePressureMatrix(const gsFunction<T> & pressFun, const gsFunctionSet<T> & deformed)
2166{
2167 try
2168 {
2169 this->_assemblePressure<true>(pressFun,deformed);
2170 m_status = ThinShellAssemblerStatus::Success;
2171 }
2172 catch (...)
2173 {
2174 m_status = ThinShellAssemblerStatus::AssemblyError;
2175 }
2176 return m_status;
2177}
2178
2179template<short_t d, class T, bool bending>
2180ThinShellAssemblerStatus gsThinShellAssembler<d, T, bending>::assemblePressureMatrix(const T pressure, const gsFunctionSet<T> & deformed)
2181{
2182 gsConstantFunction<T> pressFun(pressure,d);
2183 try
2184 {
2185 this->_assemblePressure<true>(pressFun, deformed);
2186 m_status = ThinShellAssemblerStatus::Success;
2187 }
2188 catch (...)
2189 {
2190 m_status = ThinShellAssemblerStatus::AssemblyError;
2191 }
2192 return m_status;
2193}
2194
2195template<short_t d, class T, bool bending>
2196ThinShellAssemblerStatus gsThinShellAssembler<d, T, bending>::assemblePressureVector(const gsFunction<T> & pressFun )
2197{
2198 try
2199 {
2200 this->_assemblePressure<false>(pressFun);
2201 m_status = ThinShellAssemblerStatus::Success;
2202 }
2203 catch (...)
2204 {
2205 m_status = ThinShellAssemblerStatus::AssemblyError;
2206 }
2207 return m_status;
2208}
2209
2210template<short_t d, class T, bool bending>
2211ThinShellAssemblerStatus gsThinShellAssembler<d, T, bending>::assemblePressureVector(const T pressure)
2212{
2213 gsConstantFunction<T> pressFun(pressure,d);
2214 try
2215 {
2216 this->_assemblePressure<false>(pressFun);
2217 m_status = ThinShellAssemblerStatus::Success;
2218 }
2219 catch (...)
2220 {
2221 m_status = ThinShellAssemblerStatus::AssemblyError;
2222 }
2223 return m_status;
2224}
2225
2226template<short_t d, class T, bool bending>
2227ThinShellAssemblerStatus gsThinShellAssembler<d, T, bending>::assemblePressureVector(const gsFunction<T> & pressFun, const gsFunctionSet<T> & deformed)
2228{
2229 try
2230 {
2231 this->_assemblePressure<false>(pressFun,deformed);
2232 m_status = ThinShellAssemblerStatus::Success;
2233 }
2234 catch (...)
2235 {
2236 m_status = ThinShellAssemblerStatus::AssemblyError;
2237 }
2238 return m_status;
2239}
2240
2241template<short_t d, class T, bool bending>
2242ThinShellAssemblerStatus gsThinShellAssembler<d, T, bending>::assemblePressureVector(const T pressure, const gsFunctionSet<T> & deformed)
2243{
2244 gsConstantFunction<T> pressFun(pressure,d);
2245 try
2246 {
2247 this->_assemblePressure<false>(pressFun, deformed);
2248 m_status = ThinShellAssemblerStatus::Success;
2249 }
2250 catch (...)
2251 {
2252 m_status = ThinShellAssemblerStatus::AssemblyError;
2253 }
2254 return m_status;
2255}
2256
2257template<short_t d, class T, bool bending>
2258ThinShellAssemblerStatus gsThinShellAssembler<d, T, bending>::assemblePressureVector(const gsFunctionSet<T> & deformed)
2259{
2260 try
2261 {
2262 this->assemblePressureVector(*m_pressFun,deformed);
2263 m_status = ThinShellAssemblerStatus::Success;
2264 }
2265 catch (...)
2266 {
2267 m_status = ThinShellAssemblerStatus::AssemblyError;
2268 }
2269 return m_status;
2270}
2271
2272template<short_t d, class T, bool bending>
2273ThinShellAssemblerStatus gsThinShellAssembler<d, T, bending>::assembleFoundationVector(const gsFunctionSet<T> & deformed)
2274{
2275 try
2276 {
2277 this->assembleFoundationVector(*m_foundFun,deformed);
2278 m_status = ThinShellAssemblerStatus::Success;
2279 }
2280 catch (...)
2281 {
2282 m_status = ThinShellAssemblerStatus::AssemblyError;
2283 }
2284 return m_status;
2285}
2286
2287template<short_t d, class T, bool bending>
2288ThinShellAssemblerStatus gsThinShellAssembler<d, T, bending>::assembleFoundationVector(const gsFunction<T> & foundFun, const gsFunctionSet<T> & deformed)
2289{
2290 try
2291 {
2292 this->_assembleFoundation<false>(foundFun,deformed);
2293 m_status = ThinShellAssemblerStatus::Success;
2294 }
2295 catch (...)
2296 {
2297 m_status = ThinShellAssemblerStatus::AssemblyError;
2298 }
2299 return m_status;
2300}
2301
2302template<short_t d, class T, bool bending>
2303gsMatrix<T> gsThinShellAssembler<d, T, bending>::boundaryForce(const gsFunctionSet<T> & deformed, const std::vector<patchSide> & patchSides) const
2304{
2305 return boundaryForce_impl<d, bending>(deformed,patchSides);
2306}
2307
2308template <short_t d, typename T, bool bending>
2309template <short_t _d, bool _bending>
2310typename std::enable_if<(_d==3) && _bending, gsMatrix<T> >::type
2311gsThinShellAssembler<d, T, bending>::boundaryForce_impl(const gsFunctionSet<T> & deformed, const std::vector<patchSide> & patchSides) const
2312{
2313 gsExprAssembler<T> assembler;
2314 assembler.setIntegrationElements(m_basis);
2315 space u = assembler.getSpace(*m_spaceBasis, d, 0); // last argument is the space ID
2316
2317 gsBoundaryConditions<T> bc;
2318 u.setup(bc, dirichlet::l2Projection, m_continuity);
2319
2320 gsVector<T> F(d);
2321 F.setZero();
2322 if (const gsMultiBasis<T> * mbasis = dynamic_cast<const gsMultiBasis<T>*>(&u.source()))
2323 {
2324 // Collect indices of the functions on the selected boundaries
2325 std::vector<std::unordered_set<index_t>> indices(d);
2326 gsMatrix<index_t> boundary;
2327 for (std::vector<patchSide>::const_iterator bdr = patchSides.begin(); bdr != patchSides.end(); bdr++)
2328 {
2329 boundary = mbasis->basis(bdr->patch).boundary(bdr->side());
2330 for (index_t k=0; k!=boundary.rows(); k++)
2331 for (index_t c=0; c!=d; c++)
2332 indices[c].insert(u.mapper().index(boundary.at(k),bdr->patch,c));
2333 }
2334
2335 assembler.initSystem();
2336
2337 geometryMap m_ori = assembler.getMap(m_patches);
2338 geometryMap m_def = assembler.getMap(deformed);
2339
2340 // Initialize vector
2341 // m_assembler.initVector(1);
2342
2343 gsMaterialMatrixIntegrate<T,MaterialOutput::VectorN> m_S0(m_materialMatrices,&m_patches,&deformed);
2344 gsMaterialMatrixIntegrate<T,MaterialOutput::VectorM> m_S1(m_materialMatrices,&m_patches,&deformed);
2345 auto S0 = assembler.getCoeff(m_S0);
2346 auto S1 = assembler.getCoeff(m_S1);
2347
2348 gsFunctionExpr<> mult2t("1","0","0","0","1","0","0","0","2",2);
2349 auto m_m2 = assembler.getCoeff(mult2t);
2350
2351 // this->homogenizeDirichlet();
2352
2353 auto m_N = S0.tr();
2354 auto m_Em_der = flat( jac(m_def).tr() * jac(u) ) ;
2355
2356 auto m_M = S1.tr(); // output is a column
2357 auto m_Ef_der = -( deriv2(u,sn(m_def).normalized().tr() ) + deriv2(m_def,var1(u,m_def) ) ) * reshape(m_m2,3,3); //[checked]
2358
2359 // Assemble vector (slow?)
2360 try
2361 {
2362 assembler.assemble(
2363 - ( ( m_N * m_Em_der.tr() + m_M * m_Ef_der.tr() ) * meas(m_ori) ).tr()
2364 );
2365 }
2366 catch (...)
2367 {
2368 GISMO_ERROR("Assembly of the force vector failed.");
2369 }
2370
2371 // Grab and sum control point forces on boundary indices
2372 for (index_t c = 0; c != d; c++)
2373 for (std::unordered_set<index_t>::const_iterator it = indices[c].begin(); it!=indices[c].end(); it++)
2374 F[c] += assembler.rhs().at(*it);
2375 }
2376 else
2377 GISMO_ERROR("The basis is not a gsMultiBasis!");
2378
2379 return F;
2380}
2381
2382template <short_t d, typename T, bool bending>
2383template <short_t _d, bool _bending>
2384typename std::enable_if<!(_d==3 && _bending), gsMatrix<T> >::type
2385gsThinShellAssembler<d, T, bending>::boundaryForce_impl(const gsFunctionSet<T> & deformed, const std::vector<patchSide> & patchSides) const
2386{
2387 gsExprAssembler<T> assembler;
2388 assembler.setIntegrationElements(m_basis);
2389 space u = assembler.getSpace(*m_spaceBasis, d, 0); // last argument is the space ID
2390
2391 gsBoundaryConditions<T> bc;
2392 u.setup(bc, dirichlet::l2Projection, m_continuity);
2393
2394 gsVector<T> F(d);
2395 F.setZero();
2396 if (const gsMultiBasis<T> * mbasis = dynamic_cast<const gsMultiBasis<T>*>(&u.source()))
2397 {
2398 // Collect indices of the functions on the selected boundaries
2399 std::vector<std::unordered_set<index_t>> indices(d);
2400 gsMatrix<index_t> boundary;
2401 for (std::vector<patchSide>::const_iterator bdr = patchSides.begin(); bdr != patchSides.end(); bdr++)
2402 {
2403 boundary = mbasis->basis(bdr->patch).boundary(bdr->side());
2404 for (index_t k=0; k!=boundary.rows(); k++)
2405 for (index_t c=0; c!=d; c++)
2406 indices[c].insert(u.mapper().index(boundary.at(k),bdr->patch,c));
2407 }
2408
2409 assembler.initSystem();
2410
2411 geometryMap m_ori = assembler.getMap(m_patches);
2412 geometryMap m_def = assembler.getMap(deformed);
2413
2414 // Initialize vector
2415 // m_assembler.initVector(1);
2416
2417 gsMaterialMatrixIntegrate<T,MaterialOutput::VectorN> m_S0(m_materialMatrices,&m_patches,&deformed);
2418 auto S0 = assembler.getCoeff(m_S0);
2419
2420 // this->homogenizeDirichlet();
2421
2422 auto m_N = S0.tr();
2423 auto m_Em_der = flat( jac(m_def).tr() * jac(u) ) ;
2424
2425 // Assemble vector (slow?)
2426 try
2427 {
2428 assembler.assemble(
2429 - ( ( m_N * m_Em_der.tr() ) * meas(m_ori) ).tr()
2430 );
2431 }
2432 catch (...)
2433 {
2434 GISMO_ERROR("Assembly of the force vector failed.");
2435 }
2436
2437 // Grab and sum control point forces on boundary indices
2438 for (index_t c = 0; c != d; c++)
2439 for (std::unordered_set<index_t>::const_iterator it = indices[c].begin(); it!=indices[c].end(); it++)
2440 F[c] += assembler.rhs().at(*it);
2441 }
2442 else
2443 GISMO_ERROR("The basis is not a gsMultiBasis!");
2444
2445 return F;
2446}
2447
2448template<short_t d, class T, bool bending>
2449ThinShellAssemblerStatus gsThinShellAssembler<d, T, bending>::assembleVector(const gsMatrix<T> & solVector, const bool homogenize)
2450{
2451 gsMultiPatch<T> def;
2452 constructSolution(solVector, def);
2453 return assembleVector(def,homogenize);
2454}
2455
2456template<short_t d, class T, bool bending>
2457ThinShellAssemblerStatus gsThinShellAssembler<d, T, bending>::assemble(const gsFunctionSet<T> & deformed,
2458 const bool Matrix, const bool homogenize)
2459{
2460 ThinShellAssemblerStatus status;
2461 if (Matrix)
2462 {
2463 status = assembleMatrix(deformed);
2464 if (status!=ThinShellAssemblerStatus::Success)
2465 return status;
2466 }
2467
2468 return assembleVector(deformed,homogenize);
2469}
2470template<short_t d, class T, bool bending>
2471ThinShellAssemblerStatus gsThinShellAssembler<d, T, bending>::assemble(const gsMatrix<T> & solVector,
2472 const bool Matrix, const bool homogenize)
2473{
2474 gsMultiPatch<T> def;
2475 constructSolution(solVector, def);
2476 return assemble(def,Matrix,homogenize);
2477}
2478
2479template <short_t d, class T, bool bending>
2480gsMultiPatch<T> gsThinShellAssembler<d, T, bending>::_constructSolution(const gsMatrix<T> & solVector, const gsMultiPatch<T> & /*undeformed*/) const
2481{
2482 gsMultiPatch<T> mp = m_patches;
2483 gsMultiPatch<T> displacement = constructDisplacement(solVector);
2484 for ( size_t k =0; k!=displacement.nPatches(); ++k) // Deform the geometry
2485 mp.patch(k).coefs() += displacement.patch(k).coefs();; // defG points to mp_def, therefore updated
2486
2487 return mp;
2488}
2489
2490template <short_t d, class T, bool bending>
2491gsMultiPatch<T> gsThinShellAssembler<d, T, bending>::constructSolution(const gsMatrix<T> & solVector) const
2492{
2493 return _constructSolution(solVector,m_patches);
2494}
2495
2496template <short_t d, class T, bool bending>
2497void gsThinShellAssembler<d, T, bending>::constructSolution(const gsMatrix<T> & solVector, gsMultiPatch<T> & deformed) const
2498{
2499 deformed = _constructSolution(solVector,m_patches);
2500}
2501
2502template <short_t d, class T, bool bending>
2503void gsThinShellAssembler<d, T, bending>::updateMultiPatch(const gsMatrix<T> & solVector, gsMultiPatch<T> & mp) const
2504{
2505 mp = _constructSolution(solVector,mp);
2506}
2507
2508template<short_t d, class T, bool bending>
2509T gsThinShellAssembler<d, T, bending>::getArea(const gsFunctionSet<T> & geometry)
2510{
2511 m_assembler.cleanUp();
2512 GISMO_ENSURE(m_options.hasGroup("ExprAssembler"),"The option list does not contain options with the label 'ExprAssembler'!");
2513 m_assembler.setOptions(m_options.getGroup("ExprAssembler"));
2514
2515 geometryMap G = m_assembler.getMap(geometry);
2516
2517 gsExprEvaluator<T> evaluator(m_assembler);
2518 T result = evaluator.integral(meas(G));
2519 return result;
2520}
2521
2522template<short_t d, class T, bool bending>
2523T gsThinShellAssembler<d, T, bending>::getDisplacementNorm(const gsFunctionSet<T> & deformed)
2524{
2525 m_assembler.cleanUp();
2526 GISMO_ENSURE(m_options.hasGroup("ExprAssembler"),"The option list does not contain options with the label 'ExprAssembler'!");
2527 m_assembler.setOptions(m_options.getGroup("ExprAssembler"));
2528
2529 geometryMap m_ori = m_assembler.getMap(m_patches);
2530 geometryMap m_def = m_assembler.getMap(deformed);
2531
2532 auto u = m_def - m_ori;
2533
2534 gsExprEvaluator<T> evaluator(m_assembler);
2535 T result = evaluator.integral( u.tr() * u * meas(m_def));
2536 T area = evaluator.integral(meas(m_ori));
2537
2538 return std::pow(result/area,0.5);
2539}
2540
2541template<short_t d, class T, bool bending>
2542T gsThinShellAssembler<d, T, bending>::getElasticEnergy(const gsFunctionSet<T> & deformed)
2543{
2544 m_assembler.cleanUp();
2545 GISMO_ENSURE(m_options.hasGroup("ExprAssembler"),"The option list does not contain options with the label 'ExprAssembler'!");
2546 m_assembler.setOptions(m_options.getGroup("ExprAssembler"));
2547
2548 geometryMap m_ori = m_assembler.getMap(m_patches);
2549 geometryMap m_def = m_assembler.getMap(deformed);
2550
2551 gsMaterialMatrixIntegrate<T,MaterialOutput::VectorN> m_S0(m_materialMatrices,&deformed);
2552 gsMaterialMatrixIntegrate<T,MaterialOutput::VectorM> m_S1(m_materialMatrices,&deformed);
2553 auto S0 = m_assembler.getCoeff(m_S0);
2554 auto S1 = m_assembler.getCoeff(m_S1);
2555 auto u = m_def - m_ori;
2556
2557 auto m_N = S0.tr();
2558 auto m_M = S1.tr(); // output is a column
2559
2560 gsExprEvaluator<T> evaluator(m_assembler);
2561 T result = evaluator.integral(0.5 * ( u.tr() * ( m_N + m_M ).tr() ) * meas(m_def));
2562 return result;
2563}
2564
2565template<short_t d, class T, bool bending>
2566void gsThinShellAssembler<d, T, bending>::plotSolution(std::string string, const gsMatrix<T> & solVector)
2567{
2568 m_solvector = solVector;
2569 space m_space = m_assembler.trialSpace(0);
2570 solution m_solution = m_assembler.getSolution(m_space, m_solvector);
2571 geometryMap G = m_assembler.getMap(m_patches);
2572 gsExprEvaluator<T> ev(m_assembler);
2573 ev.options().setSwitch("plot.elements", false);
2574 ev.options().setInt ("plot.npts" , 500);
2575 ev.writeParaview( m_solution, G, string);
2576}
2577
2578template<short_t d, class T, bool bending>
2579T gsThinShellAssembler<d, T, bending>::interfaceErrorC0(const gsFunctionSet<T> & deformed, const ifContainer & iFaces)
2580{
2581 geometryMap G = m_assembler.getMap(deformed);
2582 gsExprEvaluator<T> ev(m_assembler);
2583 ev.integralInterface( ( G.left() - G.right() ).sqNorm() , iFaces);
2584 ev.calcSqrt();
2585 return ev.value();
2586}
2587
2588template<short_t d, class T, bool bending>
2589T gsThinShellAssembler<d, T, bending>::interfaceErrorG1(const gsFunctionSet<T> & deformed, const ifContainer & iFaces)
2590{
2591 geometryMap G = m_assembler.getMap(deformed);
2592 gsExprEvaluator<T> ev(m_assembler);
2593 ev.integralInterface( (sn(G.left()).normalized()-sn(G.right()).normalized()).sqNorm() , iFaces);
2594 ev.calcSqrt();
2595 return ev.value();
2596}
2597
2598template<short_t d, class T, bool bending>
2599T gsThinShellAssembler<d, T, bending>::interfaceErrorNormal(const gsFunctionSet<T> & deformed, const ifContainer & iFaces)
2600{
2601 geometryMap G = m_assembler.getMap(deformed);
2602 gsExprEvaluator<T> ev(m_assembler);
2603 ev.maxInterface( (sn(G.left())-sn(G.right())).norm() , iFaces);
2604 return ev.value();
2605}
2606
2607template<short_t d, class T, bool bending>
2608T gsThinShellAssembler<d, T, bending>::interfaceErrorGaussCurvature(const gsFunctionSet<T> & deformed, const ifContainer & iFaces)
2609{
2610 geometryMap G = m_assembler.getMap(deformed);
2611 gsExprEvaluator<T> ev(m_assembler);
2612 ev.maxInterface( abs( (fform(G.left() ).inv()*fform2nd(G.left() )).det() -
2613 (fform(G.right()).inv()*fform2nd(G.right())).det() ) , iFaces);
2614 return ev.value();
2615}
2616template<short_t d, class T, bool bending>
2617T gsThinShellAssembler<d, T, bending>::interfaceErrorMeanCurvature(const gsFunctionSet<T> & deformed, const ifContainer & iFaces)
2618{
2619 geometryMap G = m_assembler.getMap(deformed);
2620 gsExprEvaluator<T> ev(m_assembler);
2621 ev.maxInterface( abs( (fform(G.left() ).inv()*fform2nd(G.left() )).trace().val() -
2622 (fform(G.right()).inv()*fform2nd(G.right())).trace().val() ) , iFaces);
2623 return ev.value();
2624}
2625template<short_t d, class T, bool bending>
2626gsMultiPatch<T> gsThinShellAssembler<d, T, bending>::constructMultiPatch(const gsMatrix<T> & solVector) const
2627{
2628 m_solvector = solVector;
2629 space m_space = m_assembler.trialSpace(0);
2630 m_space.setup(m_bcs, dirichlet::l2Projection, m_continuity);
2631 const_cast<expr::gsFeSpace<T> & >(m_space).fixedPart() = m_ddofs; //CHECK FIXEDPART
2632
2633 if (const gsMappedBasis<2,T> * mbasis = dynamic_cast<const gsMappedBasis<2,T> * >(m_spaceBasis))
2634 {
2635 gsMatrix<T> tmp;
2636 const index_t dim = m_space.dim();
2637 GISMO_ASSERT(static_cast<size_t>(dim*mbasis->size())==m_mapper.mapSize(),"Something is wrong in the sizes, basis size = "<<mbasis->size()<<" mapper size = "<<m_mapper.mapSize());
2638
2639 gsMatrix<T> cc(mbasis->size(),d);
2640 cc.setZero();
2641
2642
2643 for ( index_t p =0; p!=m_patches.nPieces(); ++p) // Deform the geometry
2644 {
2645 for (index_t c = 0; c!=dim; c++) // for all components
2646 {
2647 // loop over all basis functions (even the eliminated ones)
2648 for (size_t i = 0; i < m_mapper.patchSize(p,c); ++i)
2649 {
2650 const index_t ii = m_mapper.index(i, p, c);
2651 if ( m_mapper.is_free_index(ii) ) // DoF value is in the solVector
2652 cc(i,c) = m_solvector.at(ii);
2653 else // eliminated DoF: fill with Dirichlet data
2654 {
2655 cc(i,c) = m_ddofs.at( m_mapper.global_to_bindex(ii) );
2656 }
2657 }
2658
2659 }
2660
2661 }
2662 mbasis->global_coef_to_local_coef(cc,tmp);
2663 return mbasis->exportToPatches(tmp);
2664 }
2665 else
2666 {
2667 gsMultiPatch<T> result;
2668 // Solution vector and solution variable
2669 solution m_solution = m_assembler.getSolution(m_space, m_solvector);
2670
2671 gsMatrix<T> cc;
2672 for ( index_t p =0; p!=m_patches.nPieces(); ++p) // Deform the geometry
2673 {
2674 m_solution.extract(cc, p);
2675 result.addPatch(m_basis.basis(p).makeGeometry( give(cc) )); // defG points to mp_def, therefore updated
2676 }
2677 return result;
2678 }
2679}
2680
2681template<short_t d, class T, bool bending>
2682gsMultiPatch<T> gsThinShellAssembler<d, T, bending>::constructDisplacement(const gsMatrix<T> & solVector) const
2683{
2684 return constructMultiPatch(solVector);
2685}
2686
2687template<short_t d, class T, bool bending>
2688gsMatrix<T> gsThinShellAssembler<d, T, bending>::fullSolutionVector(const gsMatrix<T> & vector) const
2689{
2690 gsMatrix<T> solVector = vector;
2691 space m_space = m_assembler.trialSpace(0);
2692 m_space.setup(m_bcs, dirichlet::l2Projection, m_continuity);
2693 solution m_solution = m_assembler.getSolution(m_space, solVector);
2694 gsMatrix<T> result;
2695 m_solution.extractFull(result);
2696 return result.col(0);
2697}
2698
2699template<short_t d, class T, bool bending>
2700gsVector<T> gsThinShellAssembler<d, T, bending>::constructSolutionVector(const gsMultiPatch<T> & displacements) const
2701{
2702 gsVector<T> result(m_mapper.freeSize());
2703
2704 for (size_t p=0; p!=displacements.nPatches(); p++)
2705 {
2706 for (size_t dim = 0; dim!=d; dim++)
2707 {
2708 for (size_t k=0; k!=m_mapper.patchSize(p,dim); k++)
2709 {
2710 if (m_mapper.is_free(k,p,dim))
2711 {
2712 result.at(m_mapper.index(k,p,dim)) = displacements.patch(p).coefs()(k,dim);
2713 }
2714 }
2715 }
2716
2717 }
2718 return result;
2719}
2720
2721template<short_t d, class T, bool bending>
2722void gsThinShellAssembler<d, T, bending>::constructDisplacement(const gsMatrix<T> & solVector, gsMultiPatch<T> & deformed) const
2723{
2724 deformed = constructDisplacement(solVector);
2725}
2726
2727// template<short_t d, class T, bool bending>
2728// void gsThinShellAssembler<d, T, bending>::constructStresses(const gsMultiPatch<T> & deformed,
2729// gsPiecewiseFunction<T> & result,
2730// stress_type::type type) const
2731// {
2732// deformed = constructDisplacement(solVector);
2733// }
2734
2735template<short_t d, class T, bool bending>
2736gsMatrix<T> gsThinShellAssembler<d, T, bending>::computePrincipalStretches(const gsMatrix<T> & u, const gsFunctionSet<T> & deformed, const T z)
2737{
2738 // gsDebug<<"Warning: Principle Stretch computation of gsThinShellAssembler is depreciated...\n";
2739 gsMatrix<T> Z(1,1);
2740 Z.setZero();
2741 gsMatrix<T> result(3,u.cols());
2742 result.setZero();
2743 gsMatrix<T> zmat(1,1);
2744 zmat<<z;
2745 this->_getOptions();
2746
2747 m_assembler.cleanUp();
2748 GISMO_ENSURE(m_options.hasGroup("ExprAssembler"),"The option list does not contain options with the label 'ExprAssembler'!");
2749 m_assembler.setOptions(m_options.getGroup("ExprAssembler"));
2750
2751 // geometryMap m_ori = m_assembler.getMap(m_patches);
2752 // geometryMap m_def = m_assembler.getMap(*m_defpatches);
2753 // m_assembler.initSystem();
2754
2755 gsMaterialMatrixEval<T,MaterialOutput::Stretch> m_mm(m_materialMatrices,&deformed,zmat);
2756 auto mm0 = m_assembler.getCoeff(m_mm);
2757
2758 gsExprEvaluator<T> evaluator(m_assembler);
2759
2760 for (index_t k = 0; k != u.cols(); ++k)
2761 result.col(k) = evaluator.eval(mm0,u.col(k));
2762 return result;
2763}
2764
2765template <short_t d, class T, bool bending>
2766gsMatrix<T> gsThinShellAssembler<d, T, bending>::computePrincipalStresses(const gsMatrix<T> & u, const gsFunctionSet<T> & deformed, const T z)
2767{
2768 // gsDebug<<"Warning: Principle Stretch computation of gsThinShellAssembler is depreciated...\n";
2769 gsMatrix<T> Z(1,1);
2770 Z.setZero();
2771 gsMatrix<T> result(2,u.cols());
2772 result.setZero();
2773 gsMatrix<T> zmat(1,1);
2774 zmat<<z;
2775 this->_getOptions();
2776
2777 m_assembler.cleanUp();
2778 GISMO_ENSURE(m_options.hasGroup("ExprAssembler"),"The option list does not contain options with the label 'ExprAssembler'!");
2779 m_assembler.setOptions(m_options.getGroup("ExprAssembler"));
2780
2781 // geometryMap m_ori = m_assembler.getMap(m_patches);
2782 // geometryMap m_def = m_assembler.getMap(*m_defpatches);
2783 // m_assembler.initSystem();
2784
2785 gsMaterialMatrixEval<T,MaterialOutput::PStress> m_mm(m_materialMatrices,&deformed,zmat);
2786 auto mm0 = m_assembler.getCoeff(m_mm);
2787
2788 gsExprEvaluator<T> evaluator(m_assembler);
2789
2790 for (index_t k = 0; k != u.cols(); ++k)
2791 {
2792 result.col(k) = evaluator.eval(mm0,u.col(k));
2793 }
2794 return result;
2795}
2796
2797template<short_t d, class T, bool bending>
2798void gsThinShellAssembler<d, T, bending>::constructStress(const gsFunctionSet<T> & deformed,
2799 gsPiecewiseFunction<T> & result,
2800 stress_type::type type)
2801{
2802 constructStress(m_patches,deformed,result,type);
2803}
2804
2805template<short_t d, class T, bool bending>
2806void gsThinShellAssembler<d, T, bending>::constructStress(
2807 const gsFunctionSet<T> & original,
2808 const gsFunctionSet<T> & deformed,
2809 gsPiecewiseFunction<T> & result,
2810 stress_type::type type)
2811{
2812 result.clear();
2813
2814 for (size_t p = 0; p < m_patches.nPatches(); ++p )
2815 result.addPiecePointer(new gsShellStressFunction<T>(original,deformed,m_materialMatrices,p,type));
2816
2817}
2818
2819// template<short_t d, class T, bool bending>
2820// gsField<T> gsThinShellAssembler<d, T, bending>::constructStress(const gsMultiPatch<T> & deformed,
2821// stress_type::type type)
2822// {
2823// gsPiecewiseFunction<T> result;
2824// result.clear();
2825
2826// for (size_t p = 0; p < m_patches.nPatches(); ++p )
2827// result.addPiecePointer(new gsShellStressFunction<d, T, bending>(m_patches,deformed,m_materialMatrices,p,type,m_assembler));
2828
2829// gsField<T> stressField(m_patches,result, true);
2830// return stressField;
2831
2832// }
2833
2834template<short_t d, class T, bool bending>
2835void gsThinShellAssembler<d, T, bending>::projectL2_into(const gsFunction<T> & /*fun*/, gsMatrix<T>& /*result*/)
2836{
2837 // /// todo: make a projection with BCs?
2838 // /// todo: test
2839 // // this->_getOptions();
2840
2841 // m_assembler.cleanUp();
2842 // GISMO_ENSURE(m_options.hasGroup("ExprAssembler"),"The option list does not contain options with the label 'ExprAssembler'!");
2843 // m_assembler.setOptions(m_options.getGroup("ExprAssembler"));
2844
2845 // geometryMap m_ori = m_assembler.getMap(m_patches);
2846
2847 // // Initialize stystem
2848 // m_assembler.initSystem();
2849
2850 // space m_space = m_assembler.trialSpace(0);
2851 // auto function = m_assembler.getCoeff(fun, m_ori);
2852 // // auto function = m_assembler.getCoeff(fun);
2853
2854 // // assemble system
2855 // m_assembler.assemble(m_space*m_space.tr()*meas(m_ori),m_space * function*meas(m_ori));
2856
2857 // gsSparseSolver<>::uPtr solver = gsSparseSolver<T>::get( m_options.askString("Solver","CGDiagonal") );
2858 // solver->compute(m_assembler.matrix());
2859 // result = solver->solve(m_assembler.rhs());
2860}
2861
2862template<short_t d, class T, bool bending>
2863void gsThinShellAssembler<d, T, bending>::projectL2_into(const gsFunction<T> & fun, gsMultiPatch<T>& mp)
2864{
2867 gsMatrix<T> tmp = projectL2(fun);
2868 mp = m_patches;
2869
2870 // Solution vector and solution variable
2871 space m_space = m_assembler.trialSpace(0);
2872 m_space.setup(m_bcs, dirichlet::l2Projection, m_continuity);
2873 const_cast<expr::gsFeSpace<T> & >(m_space).fixedPart() = m_ddofs;
2874
2875 solution m_solution = m_assembler.getSolution(m_space, tmp);
2876
2877 gsMatrix<T> cc;
2878 for ( size_t k =0; k!=mp.nPatches(); ++k) // Deform the geometry
2879 {
2880 // // extract deformed geometry
2881 m_solution.extract(cc, k);
2882 mp.patch(k).coefs() += cc; // defG points to mp_def, therefore updated
2883 }
2884}
2885
2886
2887template<short_t d, class T, bool bending>
2888gsMatrix<T> gsThinShellAssembler<d, T, bending>::projectL2(const gsFunction<T> & fun)
2889{
2892 gsMatrix<T> result;
2893 this->projectL2_into(fun,result);
2894 return result;
2895}
2896
2897template <short_t d, class T, bool bending>
2898T gsThinShellAssembler<d, T, bending>::deformationNorm(const gsMultiPatch<T> & deformed, const gsMultiPatch<T> & original)
2899{
2900 m_assembler.cleanUp();
2901 m_assembler.setOptions(m_options);
2902
2903 geometryMap ori = m_assembler.getMap(original);
2904 geometryMap def = m_assembler.getMap(deformed);
2905
2906 gsExprEvaluator<T> evaluator(m_assembler);
2907 T result = evaluator.integral(def.sqNorm() * meas(ori));
2908 return result;
2909}
2910
2911template<short_t d, class T, bool bending>
2912void gsThinShellAssembler<d, T, bending>::_ifcTest(const T tol)
2913{
2914 geometryMap m_ori = m_assembler.getMap(m_patches);
2915 gsExprEvaluator<T> ev(m_assembler);
2916
2917 m_inPlane.clear();
2918 m_outPlane.clear();
2919
2920 for (gsBoxTopology::const_iiterator it = m_patches.topology().iBegin(); it!=m_patches.topology().iEnd(); it++)
2921 {
2922 // G1 condition
2923 ev.integralInterface( (sn(m_ori.left()).normalized()-sn(m_ori.right()).normalized()).sqNorm() );
2924 ev.calcSqrt();
2925
2926 // // Continuous normal condition
2927 // ev.maxInterface( (sn(m_ori.left())-sn(m_ori.right())).norm() );
2928 // ev.calcSqrt();
2929
2930 if (ev.value() < tol)
2931 m_inPlane.push_back(*it);
2932 else
2933 m_outPlane.push_back(*it);
2934 }
2935}
2936
2937template<short_t d, class T, bool bending>
2938bool gsThinShellAssembler<d, T, bending>::_isInPlane(const boundaryInterface & /*ifc*/, const T tol)
2939{
2940 geometryMap m_ori = m_assembler.getMap(m_patches);
2941 gsExprEvaluator<T> ev(m_assembler);
2942
2943 // G1 condition
2944 ev.integralInterface( (sn(m_ori.left()).normalized()-sn(m_ori.right()).normalized()).sqNorm() );
2945 ev.calcSqrt();
2946
2947 // // Continuous normal condition
2948 // ev.maxInterface( (sn(m_ori.left())-sn(m_ori.right())).norm() );
2949 // ev.calcSqrt();
2950
2951 return (ev.value() < tol);
2952}
2953
2954
2955}// namespace gismo
gismo::expr::gsFeSpace
Definition gsExpressions.h:973
gismo::gsBoundaryConditions
Class containing a set of boundary conditions.
Definition gsBoundaryConditions.h:342
gismo::gsConstantFunction
Class defining a globally constant function.
Definition gsConstantFunction.h:34
gismo::gsExprAssembler
Definition gsExprAssembler.h:33
gismo::gsExprEvaluator
Generic evaluator of isogeometric expressions.
Definition gsExprEvaluator.h:39
gismo::gsExprEvaluator::integral
T integral(const expr::_expr< E > &expr)
Calculates the integral of the expression expr on the whole integration domain.
Definition gsExprEvaluator.h:154
gismo::gsExprEvaluator::eval
void eval(const expr::_expr< E > &expr, gsGridIterator< T, mode, d > &git, const index_t patchInd=0)
gismo::gsExprEvaluator::writeParaview
void writeParaview(const expr::_expr< E > &expr, geometryMap G, std::string const &fn)
Creates a paraview file named fn containing valies of the.
Definition gsExprEvaluator.h:342
gismo::gsFunctionExpr
Class defining a multivariate (real or vector) function given by a string mathematical expression.
Definition gsFunctionExpr.h:52
gismo::gsFunctionSet
Interface for the set of functions defined on a domain (the total number of functions in the set equa...
Definition gsFunctionSet.h:219
gismo::gsFunctionSet::domainDim
virtual short_t domainDim() const =0
Dimension of the (source) domain.
gismo::gsFunction
A function from a n-dimensional domain to an m-dimensional image.
Definition gsFunction.h:60
gismo::gsMaterialMatrixBase
This class defines the base class for material matrices.
Definition gsMaterialMatrixBase.h:33
gismo::gsMaterialMatrixBase::uPtr
memory::unique_ptr< gsMaterialMatrixBase > uPtr
Unique pointer for gsGeometry.
Definition gsMaterialMatrixBase.h:44
gismo::gsMaterialMatrixContainer
This class serves as the evaluator of material matrices, based on gsMaterialMatrixBase.
Definition gsMaterialMatrixContainer.h:34
gismo::gsMatrix
A matrix with arbitrary coefficient type and fixed or dynamic size.
Definition gsMatrix.h:41
gismo::gsMatrix::at
T at(index_t i) const
Returns the i-th element of the vectorization of the matrix.
Definition gsMatrix.h:211
gismo::gsMultiBasis
Holds a set of patch-wise bases and their topology information.
Definition gsMultiBasis.h:37
gismo::gsMultiPatch
Container class for a set of geometry patches and their topology, that is, the interface connections ...
Definition gsMultiPatch.h:100
gismo::gsMultiPatch::patch
gsGeometry< T > & patch(size_t i) const
Return the i-th patch.
Definition gsMultiPatch.h:292
gismo::gsMultiPatch::addPatch
index_t addPatch(typename gsGeometry< T >::uPtr g)
Add a patch from a gsGeometry<T>::uPtr.
Definition gsMultiPatch.hpp:211
gismo::gsMultiPatch::nPatches
size_t nPatches() const
Number of patches.
Definition gsMultiPatch.h:274
gismo::gsOptionList
Class which holds a list of parameters/options, and provides easy access to them.
Definition gsOptionList.h:33
gismo::gsOptionList::update
void update(const gsOptionList &other, updateType type=ignoreIfUnknown)
Updates the object using the data from other.
Definition gsOptionList.cpp:253
gismo::gsOptionList::setInt
void setInt(const std::string &label, const index_t &value)
Sets an existing option label to be equal to value.
Definition gsOptionList.cpp:158
gismo::gsOptionList::wrapIntoGroup
gsOptionList wrapIntoGroup(const std::string &gn) const
Creates a new gsOptionList where all labels are wrapped into a groupname gn.
Definition gsOptionList.cpp:284
gismo::gsOptionList::setSwitch
void setSwitch(const std::string &label, const bool &value)
Sets an existing option label to be equal to value.
Definition gsOptionList.cpp:174
gismo::gsPiecewiseFunction
A function depending on an index i, typically referring to a patch/sub-domain. On each patch a differ...
Definition gsPiecewiseFunction.h:29
gismo::gsPiecewiseFunction::clear
void clear()
Clear (delete) all functions.
Definition gsPiecewiseFunction.h:148
gismo::gsShellStressFunction
Compute Cauchy stresses for a previously computed/defined displacement field. Can be pushed into gsPi...
Definition gsThinShellFunctions.h:79
gismo::gsThinShellAssembler
Assembles the system matrix and vectors for 2D and 3D shell problems, including geometric nonlinearit...
Definition gsThinShellAssembler.h:77
gismo::gsThinShellAssembler::operator=
gsThinShellAssembler & operator=(const gsThinShellAssembler &other)
Assignment operator.
Definition gsThinShellAssembler.hpp:99
gismo::gsThinShellAssembler::boundaryForce
gsMatrix< T > boundaryForce(const gsFunctionSet< T > &deformed, const std::vector< patchSide > &patchSides) const
See gsThinShellAssemblerBase for details.
Definition gsThinShellAssembler.hpp:2303
gismo::gsThinShellAssembler::assembleFoundationVector
ThinShellAssemblerStatus assembleFoundationVector(const gsFunctionSet< T > &deformed)
See gsThinShellAssemblerBase for details.
Definition gsThinShellAssembler.hpp:2273
gismo::gsThinShellAssembler::setOptions
void setOptions(gsOptionList &options)
See gsThinShellAssemblerBase for details.
Definition gsThinShellAssembler.hpp:229
gismo::gsThinShellAssembler::getDisplacementNorm
T getDisplacementNorm(const gsFunctionSet< T > &deformed)
See gsThinShellAssemblerBase for details.
Definition gsThinShellAssembler.hpp:2523
gismo::gsThinShellAssembler::constructDisplacement
gsMultiPatch< T > constructDisplacement(const gsMatrix< T > &solVector) const
See gsThinShellAssemblerBase for details.
Definition gsThinShellAssembler.hpp:2682
gismo::gsThinShellAssembler::_initialize
void _initialize()
Initializes the method.
Definition gsThinShellAssembler.hpp:243
gismo::gsThinShellAssembler::deformationNorm
T deformationNorm(const gsMultiPatch< T > &deformed, const gsMultiPatch< T > &original)
See gsThinShellAssemblerBase for details.
Definition gsThinShellAssembler.hpp:2898
gismo::gsThinShellAssembler::gsThinShellAssembler
gsThinShellAssembler()
Constructor for te shell assembler.
Definition gsThinShellAssembler.h:134
gismo::gsThinShellAssembler::computePrincipalStresses
gsMatrix< T > computePrincipalStresses(const gsMatrix< T > &points, const gsFunctionSet< T > &deformed, const T z=0)
See gsThinShellAssemblerBase for details.
Definition gsThinShellAssembler.hpp:2766
gismo::gsThinShellAssembler::assembleVector
ThinShellAssemblerStatus assembleVector(const gsFunctionSet< T > &deformed, const bool homogenize=true)
See gsThinShellAssemblerBase for details.
Definition gsThinShellAssembler.hpp:1983
gismo::gsThinShellAssembler::plotSolution
void plotSolution(std::string string, const gsMatrix< T > &solVector)
See gsThinShellAssemblerBase for details.
Definition gsThinShellAssembler.hpp:2566
gismo::gsThinShellAssembler::interfaceErrorNormal
T interfaceErrorNormal(const gsFunctionSet< T > &deformed)
See gsThinShellAssemblerBase for details.
Definition gsThinShellAssembler.h:485
gismo::gsThinShellAssembler::getArea
T getArea(const gsFunctionSet< T > &geometry)
See gsThinShellAssemblerBase for details.
Definition gsThinShellAssembler.hpp:2509
gismo::gsThinShellAssembler::interfaceErrorG1
T interfaceErrorG1(const gsFunctionSet< T > &deformed)
See gsThinShellAssemblerBase for details.
Definition gsThinShellAssembler.h:480
gismo::gsThinShellAssembler::assemblePressureVector
ThinShellAssemblerStatus assemblePressureVector(const gsFunction< T > &pressFun)
See gsThinShellAssemblerBase for details.
Definition gsThinShellAssembler.hpp:2196
gismo::gsThinShellAssembler::getElasticEnergy
T getElasticEnergy(const gsFunctionSet< T > &deformed)
See gsThinShellAssemblerBase for details.
Definition gsThinShellAssembler.hpp:2542
gismo::gsThinShellAssembler::assemble
ThinShellAssemblerStatus assemble()
See gsThinShellAssemblerBase for details.
Definition gsThinShellAssembler.hpp:1525
gismo::gsThinShellAssembler::assembleMass
ThinShellAssemblerStatus assembleMass(const bool lumped=false)
See gsThinShellAssemblerBase for details.
Definition gsThinShellAssembler.hpp:1436
gismo::gsThinShellAssembler::constructStress
void constructStress(const gsFunctionSet< T > &deformed, gsPiecewiseFunction< T > &result, stress_type::type type)
See gsThinShellAssemblerBase for details.
Definition gsThinShellAssembler.hpp:2798
gismo::gsThinShellAssembler::fullSolutionVector
gsMatrix< T > fullSolutionVector(const gsMatrix< T > &vector) const
See gsThinShellAssemblerBase for details.
Definition gsThinShellAssembler.hpp:2688
gismo::gsThinShellAssembler::assembleMatrix
ThinShellAssemblerStatus assembleMatrix(const gsFunctionSet< T > &deformed)
See gsThinShellAssemblerBase for details.
Definition gsThinShellAssembler.hpp:1706
gismo::gsThinShellAssembler::interfaceErrorGaussCurvature
T interfaceErrorGaussCurvature(const gsFunctionSet< T > &deformed)
See gsThinShellAssemblerBase for details.
Definition gsThinShellAssembler.h:490
gismo::gsThinShellAssembler::constructSolution
gsMultiPatch< T > constructSolution(const gsMatrix< T > &solVector) const
Construct deformed shell geometry from computed solution vector solVector and returns a multipatch.
Definition gsThinShellAssembler.hpp:2491
gismo::gsThinShellAssembler::projectL2_into
void projectL2_into(const gsFunction< T > &fun, gsMatrix< T > &result)
See gsThinShellAssemblerBase for details.
Definition gsThinShellAssembler.hpp:2835
gismo::gsThinShellAssembler::constructMultiPatch
gsMultiPatch< T > constructMultiPatch(const gsMatrix< T > &solVector) const
Construct solution field from computed solution vector solVector and returns a multipatch.
Definition gsThinShellAssembler.hpp:2626
gismo::gsThinShellAssembler::assembleFoundation
ThinShellAssemblerStatus assembleFoundation()
See gsThinShellAssemblerBase for details.
Definition gsThinShellAssembler.hpp:1496
gismo::gsThinShellAssembler::interfaceErrorMeanCurvature
T interfaceErrorMeanCurvature(const gsFunctionSet< T > &deformed)
See gsThinShellAssemblerBase for details.
Definition gsThinShellAssembler.h:495
gismo::gsThinShellAssembler::constructSolutionVector
gsVector< T > constructSolutionVector(const gsMultiPatch< T > &deformed) const
See gsThinShellAssemblerBase for details.
Definition gsThinShellAssembler.hpp:2700
gismo::gsThinShellAssembler::interfaceErrorC0
T interfaceErrorC0(const gsFunctionSet< T > &deformed)
See gsThinShellAssemblerBase for details.
Definition gsThinShellAssembler.h:475
gismo::gsThinShellAssembler::projectL2
gsMatrix< T > projectL2(const gsFunction< T > &fun)
See gsThinShellAssemblerBase for details.
Definition gsThinShellAssembler.hpp:2888
gismo::gsThinShellAssembler::_constructSolution
gsMultiPatch< T > _constructSolution(const gsMatrix< T > &solVector, const gsMultiPatch< T > &undeformed) const
See gsThinShellAssemblerBase for details.
Definition gsThinShellAssembler.hpp:2480
gismo::gsThinShellAssembler::homogenizeDirichlet
void homogenizeDirichlet()
See gsThinShellAssemblerBase for details.
Definition gsThinShellAssembler.hpp:1314
gismo::gsThinShellAssembler::assemblePressureMatrix
ThinShellAssemblerStatus assemblePressureMatrix(const gsFunction< T > &pressFun)
See gsThinShellAssemblerBase for details.
Definition gsThinShellAssembler.hpp:2134
gismo::gsThinShellAssembler::computePrincipalStretches
gsMatrix< T > computePrincipalStretches(const gsMatrix< T > &points, const gsFunctionSet< T > &deformed, const T z=0)
See gsThinShellAssemblerBase for details.
Definition gsThinShellAssembler.hpp:2736
gismo::gsVector
A vector with arbitrary coefficient type and fixed or dynamic size.
Definition gsVector.h:37
gismo::gsVector::at
T at(index_t i) const
Returns the i-th element of the vector.
Definition gsVector.h:177
gsBoundaryConditions.h
Provides gsBoundaryConditions class.
index_t
#define index_t
Definition gsConfig.h:32
gsConstantFunction.h
Provides declaration of ConstantFunction class.
GISMO_ERROR
#define GISMO_ERROR(message)
Definition gsDebug.h:118
GISMO_ENSURE
#define GISMO_ENSURE(cond, message)
Definition gsDebug.h:102
GISMO_ASSERT
#define GISMO_ASSERT(cond, message)
Definition gsDebug.h:89
gsFunctionExpr.h
Provides declaration of FunctionExpr class.
gsMaterialMatrixBase.h
Provides a base class for material matrices.
gsMaterialMatrixEval.h
Provides an evaluator for material matrices for thin shells.
gsMaterialMatrixIntegrate.h
Provides an evaluator for material matrices for thin shells.
gsPiecewiseFunction.h
Provides declaration of a gsPiecewiseFunction class.
gsThinShellAssembler.h
Provides linear and nonlinear assemblers for thin shells.
gismo::expr::reshape
EIGEN_STRONG_INLINE reshape_expr< E > const reshape(E const &u, index_t n, index_t m)
Reshape an expression.
Definition gsExpressions.h:1925
gismo::expr::abs
EIGEN_STRONG_INLINE abs_expr< E > abs(const E &u)
Absolute value.
Definition gsExpressions.h:4486
gismo::expr::fform2nd
EIGEN_STRONG_INLINE fform2nd_expr< T > fform2nd(const gsGeometryMap< T > &G)
The second fundamental form of G.
Definition gsExpressions.h:4523
gismo::expr::meas
EIGEN_STRONG_INLINE meas_expr< T > meas(const gsGeometryMap< T > &G)
The measure of a geometry map.
Definition gsExpressions.h:4555
gismo::expr::nv
EIGEN_STRONG_INLINE onormal_expr< T > nv(const gsGeometryMap< T > &u)
The (outer pointing) boundary normal of a geometry map.
Definition gsExpressions.h:4506
gismo
The G+Smo namespace, containing all definitions for the library.
gismo::ThinShellAssemblerStatus
ThinShellAssemblerStatus
Definition gsThinShellAssembler.h:54
gismo::ThinShellAssemblerStatus::Success
@ Success
Assembly is successful.
gismo::ThinShellAssemblerStatus::AssemblyError
@ AssemblyError
Assembly failed due to an error in the expression (e.g. overflow)
gismo::give
S give(S &x)
Definition gsMemory.h:266
gismo::boundary
Struct that defines the boundary sides and corners and types of a geometric object.
Definition gsBoundary.h:56
gismo::stress_type::type
type
Definition gsThinShellFunctions.h:39