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 GISMO_ASSERT(m_parametricForce,"The force must be defined in the parametric domain for 2D problems");
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 m_assembler.assemble(m_space * m_parforce * meas(m_ori));
1678
1679 this->_assembleWeakBCs<true>();
1680 this->_assembleWeakBCs<false>();
1681 this->_assembleWeakIfc<true>();
1682 this->_assembleWeakIfc<false>();
1683 this->_assembleNeumann();
1684
1685 // Assemble the loads
1686 if ( m_pLoads.numLoads() != 0 )
1687 {
1688 m_rhs = m_assembler.rhs();
1689 _applyLoads();
1690 }
1691
1692 m_status = ThinShellAssemblerStatus::Success;
1693 }
1694 catch (...)
1695 {
1696 m_assembler.cleanUp();
1697 m_status = ThinShellAssemblerStatus::AssemblyError;
1698 }
1699 return m_status;
1700}
1701
1702template<short_t d, class T, bool bending>
1703ThinShellAssemblerStatus gsThinShellAssembler<d, T, bending>::assembleMatrix(const gsFunctionSet<T> & deformed)
1704{
1705 return assembleMatrix_impl<d, bending>(deformed);
1706}
1707
1708template <short_t d, typename T, bool bending>
1709template <short_t _d, bool _bending>
1710typename std::enable_if<(_d==3) && _bending, ThinShellAssemblerStatus>::type
1711gsThinShellAssembler<d, T, bending>::assembleMatrix_impl(const gsFunctionSet<T> & deformed)
1712{
1713 m_assembler.cleanUp();
1714 GISMO_ENSURE(m_options.hasGroup("ExprAssembler"),"The option list does not contain options with the label 'ExprAssembler'!");
1715 m_assembler.setOptions(m_options.getGroup("ExprAssembler"));
1716
1717 geometryMap m_ori = m_assembler.getMap(m_patches);
1718 geometryMap m_def = m_assembler.getMap(deformed);
1719
1720 // Initialize matrix
1721 m_assembler.initSystem();
1722 m_assembler.initMatrix();
1723
1724 gsMaterialMatrixIntegrate<T,MaterialOutput::MatrixA> m_mmA(m_materialMatrices,&m_patches,&deformed);
1725 gsMaterialMatrixIntegrate<T,MaterialOutput::MatrixB> m_mmB(m_materialMatrices,&m_patches,&deformed);
1726 gsMaterialMatrixIntegrate<T,MaterialOutput::MatrixC> m_mmC(m_materialMatrices,&m_patches,&deformed);
1727 gsMaterialMatrixIntegrate<T,MaterialOutput::MatrixD> m_mmD(m_materialMatrices,&m_patches,&deformed);
1728 gsMaterialMatrixIntegrate<T,MaterialOutput::VectorN> m_S0(m_materialMatrices,&m_patches,&deformed);
1729 gsMaterialMatrixIntegrate<T,MaterialOutput::VectorM> m_S1(m_materialMatrices,&m_patches,&deformed);
1730 auto mmA = m_assembler.getCoeff(m_mmA);
1731 auto mmB = m_assembler.getCoeff(m_mmB);
1732 auto mmC = m_assembler.getCoeff(m_mmC);
1733 auto mmD = m_assembler.getCoeff(m_mmD);
1734 auto S0 = m_assembler.getCoeff(m_S0);
1735 auto S1 = m_assembler.getCoeff(m_S1);
1736
1737 gsFunctionExpr<T> mult2t("1","0","0","0","1","0","0","0","2",2);
1738 auto m_m2 = m_assembler.getCoeff(mult2t);
1739
1740 space m_space = m_assembler.trialSpace(0);
1741
1742 this->homogenizeDirichlet();
1743
1744 auto m_N = S0.tr();
1745 auto m_Em_der = flat( jac(m_def).tr() * jac(m_space) ) ; //[checked]
1746 auto m_Em_der2 = flatdot( jac(m_space),jac(m_space).tr(), m_N ); //[checked]
1747
1748
1749 auto m_M = S1.tr(); // output is a column
1750 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]
1751 auto m_Ef_der2 = -(flatdot2( deriv2(m_space), var1(m_space,m_def).tr(), m_M ).symmetrize()
1752 + var2deriv2(m_space,m_space,m_def, m_M ));
1753
1754 auto m_N_der = m_Em_der * reshape(mmA,3,3) + m_Ef_der * reshape(mmB,3,3);
1755 auto m_M_der = m_Em_der * reshape(mmC,3,3) + m_Ef_der * reshape(mmD,3,3);
1756
1757 try
1758 {
1759 if (m_foundInd) this->_assembleFoundation<true>(*m_foundFun,deformed);
1760 if (m_pressInd) this->_assemblePressure<true>(*m_pressFun,deformed);
1761
1762 // Assemble matrix
1763 m_assembler.assemble(
1764 (
1765 m_N_der * m_Em_der.tr()
1766 +
1767 m_Em_der2
1768 +
1769 m_M_der * m_Ef_der.tr()
1770 +
1771 m_Ef_der2
1772 ) * meas(m_ori)
1773 );
1774
1775 this->_assembleWeakBCs<true>(deformed);
1776 this->_assembleWeakIfc<true>(deformed);
1777
1778 m_status = ThinShellAssemblerStatus::Success;
1779 }
1780 catch (...)
1781 {
1782 m_assembler.cleanUp();
1783 m_status = ThinShellAssemblerStatus::AssemblyError;
1784 }
1785 return m_status;
1786}
1787
1788template <short_t d, typename T, bool bending>
1789template <short_t _d, bool _bending>
1790typename std::enable_if<!(_d==3 && _bending), ThinShellAssemblerStatus>::type
1791gsThinShellAssembler<d, T, bending>::assembleMatrix_impl(const gsFunctionSet<T> & deformed)
1792{
1793 m_assembler.cleanUp();
1794 GISMO_ENSURE(m_options.hasGroup("ExprAssembler"),"The option list does not contain options with the label 'ExprAssembler'!");
1795 m_assembler.setOptions(m_options.getGroup("ExprAssembler"));
1796
1797 geometryMap m_ori = m_assembler.getMap(m_patches);
1798 geometryMap m_def = m_assembler.getMap(deformed);
1799
1800 // Initialize matrix
1801 m_assembler.initSystem();
1802 m_assembler.initMatrix();
1803
1804 gsMaterialMatrixIntegrate<T,MaterialOutput::MatrixA> m_mmA(m_materialMatrices,&m_patches,&deformed);
1805 gsMaterialMatrixIntegrate<T,MaterialOutput::VectorN> m_S0(m_materialMatrices,&m_patches,&deformed);
1806 auto mmA = m_assembler.getCoeff(m_mmA);
1807 auto S0 = m_assembler.getCoeff(m_S0);
1808
1809 space m_space = m_assembler.trialSpace(0);
1810
1811
1812 this->homogenizeDirichlet();
1813
1814 auto m_N = S0.tr();
1815 auto m_Em_der = flat( jac(m_def).tr() * jac(m_space) ) ; //[checked]
1816 auto m_Em_der2 = flatdot( jac(m_space),jac(m_space).tr(), m_N ); //[checked]
1817
1818 auto m_N_der = m_Em_der * reshape(mmA,3,3);
1819
1820 // Assemble matrix
1821 try
1822 {
1823 if (m_foundInd) this->_assembleFoundation<true>(*m_foundFun,deformed);
1824 if (m_pressInd) this->_assemblePressure<true>(*m_pressFun,deformed);
1825
1826 m_assembler.assemble(
1827 (
1828 m_N_der * m_Em_der.tr()
1829 +
1830 m_Em_der2
1831 ) * meas(m_ori)
1832 );
1833 this->_assembleWeakBCs<true>(deformed);
1834 this->_assembleWeakIfc<true>(deformed);
1835
1836 m_status = ThinShellAssemblerStatus::Success;
1837 }
1838 catch (...) // add specific cases?
1839 {
1840 m_assembler.cleanUp();
1841 m_status = ThinShellAssemblerStatus::AssemblyError;
1842 }
1843 return m_status;
1844}
1845
1846template<short_t d, class T, bool bending>
1847ThinShellAssemblerStatus gsThinShellAssembler<d, T, bending>::assembleMatrix(const gsMatrix<T> & solVector)
1848{
1849 gsMultiPatch<T> def;
1850 constructSolution(solVector, def);
1851 return assembleMatrix(def);
1852}
1853
1854template<short_t d, class T, bool bending>
1855ThinShellAssemblerStatus gsThinShellAssembler<d, T, bending>::assembleMatrix(const gsFunctionSet<T> & deformed, const gsFunctionSet<T> & previous, gsMatrix<T> & update)
1856{
1857 return assembleMatrix_impl<d, bending>(deformed, previous, update);
1858}
1859
1860template <short_t d, typename T, bool bending>
1861template <short_t _d, bool _bending>
1862typename std::enable_if<(_d==3) && _bending, ThinShellAssemblerStatus>::type
1863gsThinShellAssembler<d, T, bending>::assembleMatrix_impl(const gsFunctionSet<T> & deformed, const gsFunctionSet<T> & previous, gsMatrix<T> & update)
1864{
1865 m_assembler.cleanUp();
1866 GISMO_ENSURE(m_options.hasGroup("ExprAssembler"),"The option list does not contain options with the label 'ExprAssembler'!");
1867 m_assembler.setOptions(m_options.getGroup("ExprAssembler"));
1868
1869 geometryMap m_ori = m_assembler.getMap(m_patches);
1870 geometryMap m_def = m_assembler.getMap(deformed);
1871 geometryMap m_prev = m_assembler.getMap(previous);
1872 // Initialize matrix
1873 m_assembler.initSystem();
1874 m_assembler.initMatrix();
1875
1876 gsMaterialMatrixIntegrate<T,MaterialOutput::MatrixA> m_mmA(m_materialMatrices,&m_patches,&deformed);
1877 gsMaterialMatrixIntegrate<T,MaterialOutput::MatrixB> m_mmB(m_materialMatrices,&m_patches,&deformed);
1878 gsMaterialMatrixIntegrate<T,MaterialOutput::MatrixC> m_mmC(m_materialMatrices,&m_patches,&deformed);
1879 gsMaterialMatrixIntegrate<T,MaterialOutput::MatrixD> m_mmD(m_materialMatrices,&m_patches,&deformed);
1880 gsMaterialMatrixIntegrate<T,MaterialOutput::VectorN> m_S0(m_materialMatrices,&m_patches,&deformed);
1881 gsMaterialMatrixIntegrate<T,MaterialOutput::VectorM> m_S1(m_materialMatrices,&m_patches,&deformed);
1882 auto mmA = m_assembler.getCoeff(m_mmA);
1883 auto mmB = m_assembler.getCoeff(m_mmB);
1884 auto mmC = m_assembler.getCoeff(m_mmC);
1885 auto mmD = m_assembler.getCoeff(m_mmD);
1886 // auto S0 = m_assembler.getCoeff(m_S0);
1887 // auto S1 = m_assembler.getCoeff(m_S1);
1888
1889 gsMaterialMatrixIntegrate<T,MaterialOutput::MatrixA> m_mmAd(m_materialMatrices,&m_patches,&previous);
1890 gsMaterialMatrixIntegrate<T,MaterialOutput::MatrixB> m_mmBd(m_materialMatrices,&m_patches,&previous);
1891 gsMaterialMatrixIntegrate<T,MaterialOutput::MatrixC> m_mmCd(m_materialMatrices,&m_patches,&previous);
1892 gsMaterialMatrixIntegrate<T,MaterialOutput::MatrixD> m_mmDd(m_materialMatrices,&m_patches,&previous);
1893 gsMaterialMatrixIntegrate<T,MaterialOutput::VectorN> m_S0d(m_materialMatrices,&m_patches,&previous);
1894 gsMaterialMatrixIntegrate<T,MaterialOutput::VectorM> m_S1d(m_materialMatrices,&m_patches,&previous);
1895 auto mmAp = m_assembler.getCoeff(m_mmAd);
1896 auto mmBp = m_assembler.getCoeff(m_mmBd);
1897 auto mmCp = m_assembler.getCoeff(m_mmCd);
1898 auto mmDp = m_assembler.getCoeff(m_mmDd);
1899 auto S0 = m_assembler.getCoeff(m_S0d);
1900 auto S1 = m_assembler.getCoeff(m_S1d);
1901
1902 gsFunctionExpr<> mult2t("1","0","0","0","1","0","0","0","2",2);
1903 auto m_m2 = m_assembler.getCoeff(mult2t);
1904
1905 space m_space = m_assembler.trialSpace(0);
1906 solution m_du = m_assembler.getSolution(m_space,update);
1907
1908 this->homogenizeDirichlet();
1909
1910 auto m_E_mc = flat( jac(m_prev).tr() * grad(m_du) ) ; //[checked]
1911 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]
1912 auto m_N_c = m_E_mc * reshape(mmAp,3,3) + m_E_fc * reshape(mmBp,3,3);
1913 auto m_M_c = m_E_mc * reshape(mmCp,3,3) + m_E_fc * reshape(mmDp,3,3);
1914
1915 auto m_N = S0.tr() + m_N_c;
1916 auto m_Em_der = flat( jac(m_def).tr() * jac(m_space) ) ; //[checked]
1917 auto m_Em_der2 = flatdot( jac(m_space),jac(m_space).tr(), m_N ); //[checked]
1918
1919 auto m_M = S1.tr() + m_M_c; // output is a column
1920 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]
1921 auto m_Ef_der2 = -(flatdot2( deriv2(m_space), var1(m_space,m_def).tr(), m_M ).symmetrize()
1922 + var2deriv2(m_space,m_space,m_def, m_M ));
1923
1924 auto m_N_der = m_Em_der * reshape(mmA,3,3) + m_Ef_der * reshape(mmB,3,3);
1925 auto m_M_der = m_Em_der * reshape(mmC,3,3) + m_Ef_der * reshape(mmD,3,3);
1926
1927 try
1928 {
1929 if (m_foundInd) this->_assembleFoundation<true>(*m_foundFun,deformed);
1930 if (m_pressInd) this->_assemblePressure<true>(*m_pressFun,deformed);
1931
1932 // Assemble matrix
1933 m_assembler.assemble(
1934 (
1935 m_N_der * m_Em_der.tr()
1936 +
1937 m_Em_der2
1938 +
1939 m_M_der * m_Ef_der.tr()
1940 +
1941 m_Ef_der2
1942 ) * meas(m_ori)
1943 );
1944 this->_assembleWeakBCs<true>(deformed);
1945 this->_assembleWeakIfc<true>(deformed);
1946
1947 m_status = ThinShellAssemblerStatus::Success;
1948 }
1949 catch (...)
1950 {
1951 m_assembler.cleanUp();
1952 m_status = ThinShellAssemblerStatus::AssemblyError;
1953 }
1954 return m_status;
1955}
1956
1957// template <short_t d, typename T, bool bending>
1958// template <short_t _d, bool _bending>
1959// typename std::enable_if<!(_d==3 && _bending), ThinShellAssemblerStatus>::type
1960// gsThinShellAssembler<d, T, bending>::assembleMatrix_impl(const gsMultiPatch<T> & deformed, const gsMultiPatch<T> & previous, gsMatrix<T> & update)
1961// {
1962// GISMO_NO_IMPLEMENTATION;
1963// }
1964
1965template<short_t d, class T, bool bending>
1966ThinShellAssemblerStatus gsThinShellAssembler<d, T, bending>::assembleMatrix(const gsMatrix<T> & solVector, const gsMatrix<T> & prevVector)
1967{
1968 // gsMultiPatch<T> deformed;
1969 // constructSolution(solVector, deformed);
1970 // assembleMatrix(deformed);
1971
1972 gsMultiPatch<T> def, it;
1973 constructSolution(solVector, def);
1974 constructSolution(prevVector, it);
1975 gsMatrix<T> update = solVector - prevVector;
1976 return assembleMatrix(def,it,update);
1977}
1978
1979template<short_t d, class T, bool bending>
1980ThinShellAssemblerStatus gsThinShellAssembler<d, T, bending>::assembleVector(const gsFunctionSet<T> & deformed, const bool homogenize)
1981{
1982 return assembleVector_impl<d, bending>(deformed,homogenize);
1983}
1984
1985template <short_t d, typename T, bool bending>
1986template <short_t _d, bool _bending>
1987typename std::enable_if<(_d==3) && _bending, ThinShellAssemblerStatus>::type
1988gsThinShellAssembler<d, T, bending>::assembleVector_impl(const gsFunctionSet<T> & deformed, const bool homogenize)
1989{
1990 m_assembler.cleanUp();
1991 GISMO_ENSURE(m_options.hasGroup("ExprAssembler"),"The option list does not contain options with the label 'ExprAssembler'!");
1992 m_assembler.setOptions(m_options.getGroup("ExprAssembler"));
1993
1994 geometryMap m_ori = m_assembler.getMap(m_patches);
1995 geometryMap m_def = m_assembler.getMap(deformed);
1996
1997 // Initialize vector
1998 m_assembler.initVector(1);
1999
2000 gsMaterialMatrixIntegrate<T,MaterialOutput::VectorN> m_S0(m_materialMatrices,&m_patches,&deformed);
2001 gsMaterialMatrixIntegrate<T,MaterialOutput::VectorM> m_S1(m_materialMatrices,&m_patches,&deformed);
2002 auto S0 = m_assembler.getCoeff(m_S0);
2003 auto S1 = m_assembler.getCoeff(m_S1);
2004
2005 gsFunctionExpr<> mult2t("1","0","0","0","1","0","0","0","2",2);
2006 auto m_m2 = m_assembler.getCoeff(mult2t);
2007
2008 space m_space = m_assembler.trialSpace(0);
2009 auto m_physforce = m_assembler.getCoeff(*m_forceFun,m_ori); // force defined in physical domain
2010 auto m_parforce = m_assembler.getCoeff(*m_forceFun); // force defined in parametric domain
2011
2012
2013 if (homogenize) this->homogenizeDirichlet();
2014 else this->_assembleDirichlet();
2015
2016 auto m_N = S0.tr();
2017 auto m_Em_der = flat( jac(m_def).tr() * jac(m_space) ) ;
2018
2019 auto m_M = S1.tr(); // output is a column
2020 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]
2021
2022 try
2023 {
2024 if (m_foundInd) this->_assembleFoundation<false>(*m_foundFun,deformed);
2025 if (m_pressInd) this->_assemblePressure<false>(*m_pressFun,deformed);
2026
2027 // Assemble vector
2029 if (m_parametricForce) // Assemble the force defined in the parameter domain
2030 m_assembler.assemble(m_space * m_parforce * meas(m_ori));
2031 else // Assemble the force defined in the physical domain
2032 m_assembler.assemble(m_space * m_physforce * meas(m_ori));
2033
2035 m_assembler.assemble(
2036 (
2037 - ( ( m_N * m_Em_der.tr() + m_M * m_Ef_der.tr() ) * meas(m_ori) ).tr()
2038 )
2039 );
2040
2041 this->_assembleWeakBCs<false>(deformed);
2042 this->_assembleWeakIfc<false>(deformed);
2043 this->_assembleNeumann();
2044
2045 // Assemble the loads
2046 if ( m_pLoads.numLoads() != 0 )
2047 {
2048 m_rhs = m_assembler.rhs();
2049 _applyLoads();
2050 }
2051
2052 m_status = ThinShellAssemblerStatus::Success;
2053 }
2054 catch (...)
2055 {
2056 m_assembler.cleanUp();
2057 m_status = ThinShellAssemblerStatus::AssemblyError;
2058 }
2059 return m_status;
2060}
2061
2062template <short_t d, typename T, bool bending>
2063template <short_t _d, bool _bending>
2064typename std::enable_if<!(_d==3 && _bending), ThinShellAssemblerStatus>::type
2065gsThinShellAssembler<d, T, bending>::assembleVector_impl(const gsFunctionSet<T> & deformed, const bool homogenize)
2066{
2067 m_assembler.cleanUp();
2068 GISMO_ENSURE(m_options.hasGroup("ExprAssembler"),"The option list does not contain options with the label 'ExprAssembler'!");
2069 m_assembler.setOptions(m_options.getGroup("ExprAssembler"));
2070
2071 geometryMap m_ori = m_assembler.getMap(m_patches);
2072 geometryMap m_def = m_assembler.getMap(deformed);
2073
2074 // Initialize vector
2075 m_assembler.initVector(1);
2076
2077 gsMaterialMatrixIntegrate<T,MaterialOutput::VectorN> m_S0(m_materialMatrices,&m_patches,&deformed);
2078 auto S0 = m_assembler.getCoeff(m_S0);
2079
2080 space m_space = m_assembler.trialSpace(0);
2081 auto m_physforce = m_assembler.getCoeff(*m_forceFun,m_ori); // force defined in physical domain
2082 auto m_parforce = m_assembler.getCoeff(*m_forceFun); // force defined in parametric domain
2083
2084 if (homogenize) this->homogenizeDirichlet();
2085 else this->_assembleDirichlet();
2086
2087 auto m_N = S0.tr();
2088 auto m_Em_der = flat( jac(m_def).tr() * jac(m_space) ) ;
2089
2090 try
2091 {
2092 if (m_foundInd) this->_assembleFoundation<false>(*m_foundFun,deformed);
2093 if (m_pressInd) this->_assemblePressure<false>(*m_pressFun,deformed);
2094
2095 // Assemble vector
2097 if (m_parametricForce) // Assemble the force defined in the parameter domain
2098 m_assembler.assemble(m_space * m_parforce * meas(m_ori));
2099 else // Assemble the force defined in the physical domain
2100 m_assembler.assemble(m_space * m_physforce * meas(m_ori));
2101
2103 m_assembler.assemble(
2104 (
2105 - ( ( m_N * m_Em_der.tr() ) * meas(m_ori) ).tr()
2106 )
2107 );
2108
2109 this->_assembleWeakBCs<false>(deformed);
2110 this->_assembleWeakIfc<false>(deformed);
2111 this->_assembleNeumann();
2112
2113 // Assemble the loads
2114 if ( m_pLoads.numLoads() != 0 )
2115 {
2116 m_rhs = m_assembler.rhs();
2117 _applyLoads();
2118 }
2119
2120 m_status = ThinShellAssemblerStatus::Success;
2121 }
2122 catch (...)
2123 {
2124 m_assembler.cleanUp();
2125 m_status = ThinShellAssemblerStatus::AssemblyError;
2126 }
2127 return m_status;
2128}
2129
2130template<short_t d, class T, bool bending>
2131ThinShellAssemblerStatus gsThinShellAssembler<d, T, bending>::assemblePressureMatrix(const gsFunction<T> & pressFun)
2132{
2133 try
2134 {
2135 this->_assemblePressure<true>(pressFun);
2136 m_status = ThinShellAssemblerStatus::Success;
2137 }
2138 catch (...)
2139 {
2140 m_status = ThinShellAssemblerStatus::AssemblyError;
2141 }
2142 return m_status;
2143}
2144
2145template<short_t d, class T, bool bending>
2146ThinShellAssemblerStatus gsThinShellAssembler<d, T, bending>::assemblePressureMatrix(const T pressure)
2147{
2148 gsConstantFunction<T> pressFun(pressure,d);
2149 try
2150 {
2151 this->_assemblePressure<true>(pressFun);
2152 m_status = ThinShellAssemblerStatus::Success;
2153 }
2154 catch (...)
2155 {
2156 m_status = ThinShellAssemblerStatus::AssemblyError;
2157 }
2158 return m_status;
2159}
2160
2161template<short_t d, class T, bool bending>
2162ThinShellAssemblerStatus gsThinShellAssembler<d, T, bending>::assemblePressureMatrix(const gsFunction<T> & pressFun, const gsFunctionSet<T> & deformed)
2163{
2164 try
2165 {
2166 this->_assemblePressure<true>(pressFun,deformed);
2167 m_status = ThinShellAssemblerStatus::Success;
2168 }
2169 catch (...)
2170 {
2171 m_status = ThinShellAssemblerStatus::AssemblyError;
2172 }
2173 return m_status;
2174}
2175
2176template<short_t d, class T, bool bending>
2177ThinShellAssemblerStatus gsThinShellAssembler<d, T, bending>::assemblePressureMatrix(const T pressure, const gsFunctionSet<T> & deformed)
2178{
2179 gsConstantFunction<T> pressFun(pressure,d);
2180 try
2181 {
2182 this->_assemblePressure<true>(pressFun, deformed);
2183 m_status = ThinShellAssemblerStatus::Success;
2184 }
2185 catch (...)
2186 {
2187 m_status = ThinShellAssemblerStatus::AssemblyError;
2188 }
2189 return m_status;
2190}
2191
2192template<short_t d, class T, bool bending>
2193ThinShellAssemblerStatus gsThinShellAssembler<d, T, bending>::assemblePressureVector(const gsFunction<T> & pressFun )
2194{
2195 try
2196 {
2197 this->_assemblePressure<false>(pressFun);
2198 m_status = ThinShellAssemblerStatus::Success;
2199 }
2200 catch (...)
2201 {
2202 m_status = ThinShellAssemblerStatus::AssemblyError;
2203 }
2204 return m_status;
2205}
2206
2207template<short_t d, class T, bool bending>
2208ThinShellAssemblerStatus gsThinShellAssembler<d, T, bending>::assemblePressureVector(const T pressure)
2209{
2210 gsConstantFunction<T> pressFun(pressure,d);
2211 try
2212 {
2213 this->_assemblePressure<false>(pressFun);
2214 m_status = ThinShellAssemblerStatus::Success;
2215 }
2216 catch (...)
2217 {
2218 m_status = ThinShellAssemblerStatus::AssemblyError;
2219 }
2220 return m_status;
2221}
2222
2223template<short_t d, class T, bool bending>
2224ThinShellAssemblerStatus gsThinShellAssembler<d, T, bending>::assemblePressureVector(const gsFunction<T> & pressFun, const gsFunctionSet<T> & deformed)
2225{
2226 try
2227 {
2228 this->_assemblePressure<false>(pressFun,deformed);
2229 m_status = ThinShellAssemblerStatus::Success;
2230 }
2231 catch (...)
2232 {
2233 m_status = ThinShellAssemblerStatus::AssemblyError;
2234 }
2235 return m_status;
2236}
2237
2238template<short_t d, class T, bool bending>
2239ThinShellAssemblerStatus gsThinShellAssembler<d, T, bending>::assemblePressureVector(const T pressure, const gsFunctionSet<T> & deformed)
2240{
2241 gsConstantFunction<T> pressFun(pressure,d);
2242 try
2243 {
2244 this->_assemblePressure<false>(pressFun, deformed);
2245 m_status = ThinShellAssemblerStatus::Success;
2246 }
2247 catch (...)
2248 {
2249 m_status = ThinShellAssemblerStatus::AssemblyError;
2250 }
2251 return m_status;
2252}
2253
2254template<short_t d, class T, bool bending>
2255ThinShellAssemblerStatus gsThinShellAssembler<d, T, bending>::assemblePressureVector(const gsFunctionSet<T> & deformed)
2256{
2257 try
2258 {
2259 this->assemblePressureVector(*m_pressFun,deformed);
2260 m_status = ThinShellAssemblerStatus::Success;
2261 }
2262 catch (...)
2263 {
2264 m_status = ThinShellAssemblerStatus::AssemblyError;
2265 }
2266 return m_status;
2267}
2268
2269template<short_t d, class T, bool bending>
2270ThinShellAssemblerStatus gsThinShellAssembler<d, T, bending>::assembleFoundationVector(const gsFunctionSet<T> & deformed)
2271{
2272 try
2273 {
2274 this->assembleFoundationVector(*m_foundFun,deformed);
2275 m_status = ThinShellAssemblerStatus::Success;
2276 }
2277 catch (...)
2278 {
2279 m_status = ThinShellAssemblerStatus::AssemblyError;
2280 }
2281 return m_status;
2282}
2283
2284template<short_t d, class T, bool bending>
2285ThinShellAssemblerStatus gsThinShellAssembler<d, T, bending>::assembleFoundationVector(const gsFunction<T> & foundFun, const gsFunctionSet<T> & deformed)
2286{
2287 try
2288 {
2289 this->_assembleFoundation<false>(foundFun,deformed);
2290 m_status = ThinShellAssemblerStatus::Success;
2291 }
2292 catch (...)
2293 {
2294 m_status = ThinShellAssemblerStatus::AssemblyError;
2295 }
2296 return m_status;
2297}
2298
2299template<short_t d, class T, bool bending>
2300gsMatrix<T> gsThinShellAssembler<d, T, bending>::boundaryForce(const gsFunctionSet<T> & deformed, const std::vector<patchSide> & patchSides) const
2301{
2302 return boundaryForce_impl<d, bending>(deformed,patchSides);
2303}
2304
2305template <short_t d, typename T, bool bending>
2306template <short_t _d, bool _bending>
2307typename std::enable_if<(_d==3) && _bending, gsMatrix<T> >::type
2308gsThinShellAssembler<d, T, bending>::boundaryForce_impl(const gsFunctionSet<T> & deformed, const std::vector<patchSide> & patchSides) const
2309{
2310 gsExprAssembler<T> assembler;
2311 assembler.setIntegrationElements(m_basis);
2312 space u = assembler.getSpace(*m_spaceBasis, d, 0); // last argument is the space ID
2313
2314 gsBoundaryConditions<T> bc;
2315 u.setup(bc, dirichlet::l2Projection, m_continuity);
2316
2317 gsVector<T> F(d);
2318 F.setZero();
2319 if (const gsMultiBasis<T> * mbasis = dynamic_cast<const gsMultiBasis<T>*>(&u.source()))
2320 {
2321 // Collect indices of the functions on the selected boundaries
2322 std::vector<std::unordered_set<index_t>> indices(d);
2323 gsMatrix<index_t> boundary;
2324 for (std::vector<patchSide>::const_iterator bdr = patchSides.begin(); bdr != patchSides.end(); bdr++)
2325 {
2326 boundary = mbasis->basis(bdr->patch).boundary(bdr->side());
2327 for (index_t k=0; k!=boundary.rows(); k++)
2328 for (index_t c=0; c!=d; c++)
2329 indices[c].insert(u.mapper().index(boundary.at(k),bdr->patch,c));
2330 }
2331
2332 assembler.initSystem();
2333
2334 geometryMap m_ori = assembler.getMap(m_patches);
2335 geometryMap m_def = assembler.getMap(deformed);
2336
2337 // Initialize vector
2338 // m_assembler.initVector(1);
2339
2340 gsMaterialMatrixIntegrate<T,MaterialOutput::VectorN> m_S0(m_materialMatrices,&m_patches,&deformed);
2341 gsMaterialMatrixIntegrate<T,MaterialOutput::VectorM> m_S1(m_materialMatrices,&m_patches,&deformed);
2342 auto S0 = assembler.getCoeff(m_S0);
2343 auto S1 = assembler.getCoeff(m_S1);
2344
2345 gsFunctionExpr<> mult2t("1","0","0","0","1","0","0","0","2",2);
2346 auto m_m2 = assembler.getCoeff(mult2t);
2347
2348 // this->homogenizeDirichlet();
2349
2350 auto m_N = S0.tr();
2351 auto m_Em_der = flat( jac(m_def).tr() * jac(u) ) ;
2352
2353 auto m_M = S1.tr(); // output is a column
2354 auto m_Ef_der = -( deriv2(u,sn(m_def).normalized().tr() ) + deriv2(m_def,var1(u,m_def) ) ) * reshape(m_m2,3,3); //[checked]
2355
2356 // Assemble vector (slow?)
2357 try
2358 {
2359 assembler.assemble(
2360 - ( ( m_N * m_Em_der.tr() + m_M * m_Ef_der.tr() ) * meas(m_ori) ).tr()
2361 );
2362 }
2363 catch (...)
2364 {
2365 GISMO_ERROR("Assembly of the force vector failed.");
2366 }
2367
2368 // Grab and sum control point forces on boundary indices
2369 for (index_t c = 0; c != d; c++)
2370 for (std::unordered_set<index_t>::const_iterator it = indices[c].begin(); it!=indices[c].end(); it++)
2371 F[c] += assembler.rhs().at(*it);
2372 }
2373 else
2374 GISMO_ERROR("The basis is not a gsMultiBasis!");
2375
2376 return F;
2377}
2378
2379template <short_t d, typename T, bool bending>
2380template <short_t _d, bool _bending>
2381typename std::enable_if<!(_d==3 && _bending), gsMatrix<T> >::type
2382gsThinShellAssembler<d, T, bending>::boundaryForce_impl(const gsFunctionSet<T> & deformed, const std::vector<patchSide> & patchSides) const
2383{
2384 gsExprAssembler<T> assembler;
2385 assembler.setIntegrationElements(m_basis);
2386 space u = assembler.getSpace(*m_spaceBasis, d, 0); // last argument is the space ID
2387
2388 gsBoundaryConditions<T> bc;
2389 u.setup(bc, dirichlet::l2Projection, m_continuity);
2390
2391 gsVector<T> F(d);
2392 F.setZero();
2393 if (const gsMultiBasis<T> * mbasis = dynamic_cast<const gsMultiBasis<T>*>(&u.source()))
2394 {
2395 // Collect indices of the functions on the selected boundaries
2396 std::vector<std::unordered_set<index_t>> indices(d);
2397 gsMatrix<index_t> boundary;
2398 for (std::vector<patchSide>::const_iterator bdr = patchSides.begin(); bdr != patchSides.end(); bdr++)
2399 {
2400 boundary = mbasis->basis(bdr->patch).boundary(bdr->side());
2401 for (index_t k=0; k!=boundary.rows(); k++)
2402 for (index_t c=0; c!=d; c++)
2403 indices[c].insert(u.mapper().index(boundary.at(k),bdr->patch,c));
2404 }
2405
2406 assembler.initSystem();
2407
2408 geometryMap m_ori = assembler.getMap(m_patches);
2409 geometryMap m_def = assembler.getMap(deformed);
2410
2411 // Initialize vector
2412 // m_assembler.initVector(1);
2413
2414 gsMaterialMatrixIntegrate<T,MaterialOutput::VectorN> m_S0(m_materialMatrices,&m_patches,&deformed);
2415 auto S0 = assembler.getCoeff(m_S0);
2416
2417 // this->homogenizeDirichlet();
2418
2419 auto m_N = S0.tr();
2420 auto m_Em_der = flat( jac(m_def).tr() * jac(u) ) ;
2421
2422 // Assemble vector (slow?)
2423 try
2424 {
2425 assembler.assemble(
2426 - ( ( m_N * m_Em_der.tr() ) * meas(m_ori) ).tr()
2427 );
2428 }
2429 catch (...)
2430 {
2431 GISMO_ERROR("Assembly of the force vector failed.");
2432 }
2433
2434 // Grab and sum control point forces on boundary indices
2435 for (index_t c = 0; c != d; c++)
2436 for (std::unordered_set<index_t>::const_iterator it = indices[c].begin(); it!=indices[c].end(); it++)
2437 F[c] += assembler.rhs().at(*it);
2438 }
2439 else
2440 GISMO_ERROR("The basis is not a gsMultiBasis!");
2441
2442 return F;
2443}
2444
2445template<short_t d, class T, bool bending>
2446ThinShellAssemblerStatus gsThinShellAssembler<d, T, bending>::assembleVector(const gsMatrix<T> & solVector, const bool homogenize)
2447{
2448 gsMultiPatch<T> def;
2449 constructSolution(solVector, def);
2450 return assembleVector(def,homogenize);
2451}
2452
2453template<short_t d, class T, bool bending>
2454ThinShellAssemblerStatus gsThinShellAssembler<d, T, bending>::assemble(const gsFunctionSet<T> & deformed,
2455 const bool Matrix, const bool homogenize)
2456{
2457 ThinShellAssemblerStatus status;
2458 if (Matrix)
2459 {
2460 status = assembleMatrix(deformed);
2461 if (status!=ThinShellAssemblerStatus::Success)
2462 return status;
2463 }
2464
2465 return assembleVector(deformed,homogenize);
2466}
2467template<short_t d, class T, bool bending>
2468ThinShellAssemblerStatus gsThinShellAssembler<d, T, bending>::assemble(const gsMatrix<T> & solVector,
2469 const bool Matrix, const bool homogenize)
2470{
2471 gsMultiPatch<T> def;
2472 constructSolution(solVector, def);
2473 return assemble(def,Matrix,homogenize);
2474}
2475
2476template <short_t d, class T, bool bending>
2477gsMultiPatch<T> gsThinShellAssembler<d, T, bending>::_constructSolution(const gsMatrix<T> & solVector, const gsMultiPatch<T> & /*undeformed*/) const
2478{
2479 gsMultiPatch<T> mp = m_patches;
2480 gsMultiPatch<T> displacement = constructDisplacement(solVector);
2481 for ( size_t k =0; k!=displacement.nPatches(); ++k) // Deform the geometry
2482 mp.patch(k).coefs() += displacement.patch(k).coefs();; // defG points to mp_def, therefore updated
2483
2484 return mp;
2485}
2486
2487template <short_t d, class T, bool bending>
2488gsMultiPatch<T> gsThinShellAssembler<d, T, bending>::constructSolution(const gsMatrix<T> & solVector) const
2489{
2490 return _constructSolution(solVector,m_patches);
2491}
2492
2493template <short_t d, class T, bool bending>
2494void gsThinShellAssembler<d, T, bending>::constructSolution(const gsMatrix<T> & solVector, gsMultiPatch<T> & deformed) const
2495{
2496 deformed = _constructSolution(solVector,m_patches);
2497}
2498
2499template <short_t d, class T, bool bending>
2500void gsThinShellAssembler<d, T, bending>::updateMultiPatch(const gsMatrix<T> & solVector, gsMultiPatch<T> & mp) const
2501{
2502 mp = _constructSolution(solVector,mp);
2503}
2504
2505template<short_t d, class T, bool bending>
2506T gsThinShellAssembler<d, T, bending>::getArea(const gsFunctionSet<T> & geometry)
2507{
2508 m_assembler.cleanUp();
2509 GISMO_ENSURE(m_options.hasGroup("ExprAssembler"),"The option list does not contain options with the label 'ExprAssembler'!");
2510 m_assembler.setOptions(m_options.getGroup("ExprAssembler"));
2511
2512 geometryMap G = m_assembler.getMap(geometry);
2513
2514 gsExprEvaluator<T> evaluator(m_assembler);
2515 T result = evaluator.integral(meas(G));
2516 return result;
2517}
2518
2519template<short_t d, class T, bool bending>
2520T gsThinShellAssembler<d, T, bending>::getDisplacementNorm(const gsFunctionSet<T> & deformed)
2521{
2522 m_assembler.cleanUp();
2523 GISMO_ENSURE(m_options.hasGroup("ExprAssembler"),"The option list does not contain options with the label 'ExprAssembler'!");
2524 m_assembler.setOptions(m_options.getGroup("ExprAssembler"));
2525
2526 geometryMap m_ori = m_assembler.getMap(m_patches);
2527 geometryMap m_def = m_assembler.getMap(deformed);
2528
2529 auto u = m_def - m_ori;
2530
2531 gsExprEvaluator<T> evaluator(m_assembler);
2532 T result = evaluator.integral( u.tr() * u * meas(m_def));
2533 T area = evaluator.integral(meas(m_ori));
2534
2535 return std::pow(result/area,0.5);
2536}
2537
2538template<short_t d, class T, bool bending>
2539T gsThinShellAssembler<d, T, bending>::getElasticEnergy(const gsFunctionSet<T> & deformed)
2540{
2541 m_assembler.cleanUp();
2542 GISMO_ENSURE(m_options.hasGroup("ExprAssembler"),"The option list does not contain options with the label 'ExprAssembler'!");
2543 m_assembler.setOptions(m_options.getGroup("ExprAssembler"));
2544
2545 geometryMap m_ori = m_assembler.getMap(m_patches);
2546 geometryMap m_def = m_assembler.getMap(deformed);
2547
2548 gsMaterialMatrixIntegrate<T,MaterialOutput::VectorN> m_S0(m_materialMatrices,&deformed);
2549 gsMaterialMatrixIntegrate<T,MaterialOutput::VectorM> m_S1(m_materialMatrices,&deformed);
2550 auto S0 = m_assembler.getCoeff(m_S0);
2551 auto S1 = m_assembler.getCoeff(m_S1);
2552 auto u = m_def - m_ori;
2553
2554 auto m_N = S0.tr();
2555 auto m_M = S1.tr(); // output is a column
2556
2557 gsExprEvaluator<T> evaluator(m_assembler);
2558 T result = evaluator.integral(0.5 * ( u.tr() * ( m_N + m_M ).tr() ) * meas(m_def));
2559 return result;
2560}
2561
2562template<short_t d, class T, bool bending>
2563void gsThinShellAssembler<d, T, bending>::plotSolution(std::string string, const gsMatrix<T> & solVector)
2564{
2565 m_solvector = solVector;
2566 space m_space = m_assembler.trialSpace(0);
2567 solution m_solution = m_assembler.getSolution(m_space, m_solvector);
2568 geometryMap G = m_assembler.getMap(m_patches);
2569 gsExprEvaluator<T> ev(m_assembler);
2570 ev.options().setSwitch("plot.elements", false);
2571 ev.options().setInt ("plot.npts" , 500);
2572 ev.writeParaview( m_solution, G, string);
2573}
2574
2575template<short_t d, class T, bool bending>
2576T gsThinShellAssembler<d, T, bending>::interfaceErrorC0(const gsFunctionSet<T> & deformed, const ifContainer & iFaces)
2577{
2578 geometryMap G = m_assembler.getMap(deformed);
2579 gsExprEvaluator<T> ev(m_assembler);
2580 ev.integralInterface( ( G.left() - G.right() ).sqNorm() , iFaces);
2581 ev.calcSqrt();
2582 return ev.value();
2583}
2584
2585template<short_t d, class T, bool bending>
2586T gsThinShellAssembler<d, T, bending>::interfaceErrorG1(const gsFunctionSet<T> & deformed, const ifContainer & iFaces)
2587{
2588 geometryMap G = m_assembler.getMap(deformed);
2589 gsExprEvaluator<T> ev(m_assembler);
2590 ev.integralInterface( (sn(G.left()).normalized()-sn(G.right()).normalized()).sqNorm() , iFaces);
2591 ev.calcSqrt();
2592 return ev.value();
2593}
2594
2595template<short_t d, class T, bool bending>
2596T gsThinShellAssembler<d, T, bending>::interfaceErrorNormal(const gsFunctionSet<T> & deformed, const ifContainer & iFaces)
2597{
2598 geometryMap G = m_assembler.getMap(deformed);
2599 gsExprEvaluator<T> ev(m_assembler);
2600 ev.maxInterface( (sn(G.left())-sn(G.right())).norm() , iFaces);
2601 return ev.value();
2602}
2603
2604template<short_t d, class T, bool bending>
2605T gsThinShellAssembler<d, T, bending>::interfaceErrorGaussCurvature(const gsFunctionSet<T> & deformed, const ifContainer & iFaces)
2606{
2607 geometryMap G = m_assembler.getMap(deformed);
2608 gsExprEvaluator<T> ev(m_assembler);
2609 ev.maxInterface( abs( (fform(G.left() ).inv()*fform2nd(G.left() )).det() -
2610 (fform(G.right()).inv()*fform2nd(G.right())).det() ) , iFaces);
2611 return ev.value();
2612}
2613template<short_t d, class T, bool bending>
2614T gsThinShellAssembler<d, T, bending>::interfaceErrorMeanCurvature(const gsFunctionSet<T> & deformed, const ifContainer & iFaces)
2615{
2616 geometryMap G = m_assembler.getMap(deformed);
2617 gsExprEvaluator<T> ev(m_assembler);
2618 ev.maxInterface( abs( (fform(G.left() ).inv()*fform2nd(G.left() )).trace().val() -
2619 (fform(G.right()).inv()*fform2nd(G.right())).trace().val() ) , iFaces);
2620 return ev.value();
2621}
2622template<short_t d, class T, bool bending>
2623gsMultiPatch<T> gsThinShellAssembler<d, T, bending>::constructMultiPatch(const gsMatrix<T> & solVector) const
2624{
2625 m_solvector = solVector;
2626 space m_space = m_assembler.trialSpace(0);
2627 m_space.setup(m_bcs, dirichlet::l2Projection, m_continuity);
2628 const_cast<expr::gsFeSpace<T> & >(m_space).fixedPart() = m_ddofs; //CHECK FIXEDPART
2629
2630 if (const gsMappedBasis<2,T> * mbasis = dynamic_cast<const gsMappedBasis<2,T> * >(m_spaceBasis))
2631 {
2632 gsMatrix<T> tmp;
2633 const index_t dim = m_space.dim();
2634 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());
2635
2636 gsMatrix<T> cc(mbasis->size(),d);
2637 cc.setZero();
2638
2639
2640 for ( index_t p =0; p!=m_patches.nPieces(); ++p) // Deform the geometry
2641 {
2642 for (index_t c = 0; c!=dim; c++) // for all components
2643 {
2644 // loop over all basis functions (even the eliminated ones)
2645 for (size_t i = 0; i < m_mapper.patchSize(p,c); ++i)
2646 {
2647 const index_t ii = m_mapper.index(i, p, c);
2648 if ( m_mapper.is_free_index(ii) ) // DoF value is in the solVector
2649 cc(i,c) = m_solvector.at(ii);
2650 else // eliminated DoF: fill with Dirichlet data
2651 {
2652 cc(i,c) = m_ddofs.at( m_mapper.global_to_bindex(ii) );
2653 }
2654 }
2655
2656 }
2657
2658 }
2659 mbasis->global_coef_to_local_coef(cc,tmp);
2660 return mbasis->exportToPatches(tmp);
2661 }
2662 else
2663 {
2664 gsMultiPatch<T> result;
2665 // Solution vector and solution variable
2666 solution m_solution = m_assembler.getSolution(m_space, m_solvector);
2667
2668 gsMatrix<T> cc;
2669 for ( index_t p =0; p!=m_patches.nPieces(); ++p) // Deform the geometry
2670 {
2671 m_solution.extract(cc, p);
2672 result.addPatch(m_basis.basis(p).makeGeometry( give(cc) )); // defG points to mp_def, therefore updated
2673 }
2674 return result;
2675 }
2676}
2677
2678template<short_t d, class T, bool bending>
2679gsMultiPatch<T> gsThinShellAssembler<d, T, bending>::constructDisplacement(const gsMatrix<T> & solVector) const
2680{
2681 return constructMultiPatch(solVector);
2682}
2683
2684template<short_t d, class T, bool bending>
2685gsMatrix<T> gsThinShellAssembler<d, T, bending>::fullSolutionVector(const gsMatrix<T> & vector) const
2686{
2687 gsMatrix<T> solVector = vector;
2688 space m_space = m_assembler.trialSpace(0);
2689 m_space.setup(m_bcs, dirichlet::l2Projection, m_continuity);
2690 solution m_solution = m_assembler.getSolution(m_space, solVector);
2691 gsMatrix<T> result;
2692 m_solution.extractFull(result);
2693 return result.col(0);
2694}
2695
2696template<short_t d, class T, bool bending>
2697gsVector<T> gsThinShellAssembler<d, T, bending>::constructSolutionVector(const gsMultiPatch<T> & displacements) const
2698{
2699 gsVector<T> result(m_mapper.freeSize());
2700
2701 for (size_t p=0; p!=displacements.nPatches(); p++)
2702 {
2703 for (size_t dim = 0; dim!=d; dim++)
2704 {
2705 for (size_t k=0; k!=m_mapper.patchSize(p,dim); k++)
2706 {
2707 if (m_mapper.is_free(k,p,dim))
2708 {
2709 result.at(m_mapper.index(k,p,dim)) = displacements.patch(p).coefs()(k,dim);
2710 }
2711 }
2712 }
2713
2714 }
2715 return result;
2716}
2717
2718template<short_t d, class T, bool bending>
2719void gsThinShellAssembler<d, T, bending>::constructDisplacement(const gsMatrix<T> & solVector, gsMultiPatch<T> & deformed) const
2720{
2721 deformed = constructDisplacement(solVector);
2722}
2723
2724// template<short_t d, class T, bool bending>
2725// void gsThinShellAssembler<d, T, bending>::constructStresses(const gsMultiPatch<T> & deformed,
2726// gsPiecewiseFunction<T> & result,
2727// stress_type::type type) const
2728// {
2729// deformed = constructDisplacement(solVector);
2730// }
2731
2732template<short_t d, class T, bool bending>
2733gsMatrix<T> gsThinShellAssembler<d, T, bending>::computePrincipalStretches(const gsMatrix<T> & u, const gsFunctionSet<T> & deformed, const T z)
2734{
2735 // gsDebug<<"Warning: Principle Stretch computation of gsThinShellAssembler is depreciated...\n";
2736 gsMatrix<T> Z(1,1);
2737 Z.setZero();
2738 gsMatrix<T> result(3,u.cols());
2739 result.setZero();
2740 gsMatrix<T> zmat(1,1);
2741 zmat<<z;
2742 this->_getOptions();
2743
2744 m_assembler.cleanUp();
2745 GISMO_ENSURE(m_options.hasGroup("ExprAssembler"),"The option list does not contain options with the label 'ExprAssembler'!");
2746 m_assembler.setOptions(m_options.getGroup("ExprAssembler"));
2747
2748 // geometryMap m_ori = m_assembler.getMap(m_patches);
2749 // geometryMap m_def = m_assembler.getMap(*m_defpatches);
2750 // m_assembler.initSystem();
2751
2752 gsMaterialMatrixEval<T,MaterialOutput::Stretch> m_mm(m_materialMatrices,&deformed,zmat);
2753 auto mm0 = m_assembler.getCoeff(m_mm);
2754
2755 gsExprEvaluator<T> evaluator(m_assembler);
2756
2757 for (index_t k = 0; k != u.cols(); ++k)
2758 result.col(k) = evaluator.eval(mm0,u.col(k));
2759 return result;
2760}
2761
2762template <short_t d, class T, bool bending>
2763gsMatrix<T> gsThinShellAssembler<d, T, bending>::computePrincipalStresses(const gsMatrix<T> & u, const gsFunctionSet<T> & deformed, const T z)
2764{
2765 // gsDebug<<"Warning: Principle Stretch computation of gsThinShellAssembler is depreciated...\n";
2766 gsMatrix<T> Z(1,1);
2767 Z.setZero();
2768 gsMatrix<T> result(2,u.cols());
2769 result.setZero();
2770 gsMatrix<T> zmat(1,1);
2771 zmat<<z;
2772 this->_getOptions();
2773
2774 m_assembler.cleanUp();
2775 GISMO_ENSURE(m_options.hasGroup("ExprAssembler"),"The option list does not contain options with the label 'ExprAssembler'!");
2776 m_assembler.setOptions(m_options.getGroup("ExprAssembler"));
2777
2778 // geometryMap m_ori = m_assembler.getMap(m_patches);
2779 // geometryMap m_def = m_assembler.getMap(*m_defpatches);
2780 // m_assembler.initSystem();
2781
2782 gsMaterialMatrixEval<T,MaterialOutput::PStress> m_mm(m_materialMatrices,&deformed,zmat);
2783 auto mm0 = m_assembler.getCoeff(m_mm);
2784
2785 gsExprEvaluator<T> evaluator(m_assembler);
2786
2787 for (index_t k = 0; k != u.cols(); ++k)
2788 {
2789 result.col(k) = evaluator.eval(mm0,u.col(k));
2790 }
2791 return result;
2792}
2793
2794template<short_t d, class T, bool bending>
2795void gsThinShellAssembler<d, T, bending>::constructStress(const gsFunctionSet<T> & deformed,
2796 gsPiecewiseFunction<T> & result,
2797 stress_type::type type)
2798{
2799 constructStress(m_patches,deformed,result,type);
2800}
2801
2802template<short_t d, class T, bool bending>
2803void gsThinShellAssembler<d, T, bending>::constructStress(
2804 const gsFunctionSet<T> & original,
2805 const gsFunctionSet<T> & deformed,
2806 gsPiecewiseFunction<T> & result,
2807 stress_type::type type)
2808{
2809 result.clear();
2810
2811 for (size_t p = 0; p < m_patches.nPatches(); ++p )
2812 result.addPiecePointer(new gsShellStressFunction<T>(original,deformed,m_materialMatrices,p,type));
2813
2814}
2815
2816// template<short_t d, class T, bool bending>
2817// gsField<T> gsThinShellAssembler<d, T, bending>::constructStress(const gsMultiPatch<T> & deformed,
2818// stress_type::type type)
2819// {
2820// gsPiecewiseFunction<T> result;
2821// result.clear();
2822
2823// for (size_t p = 0; p < m_patches.nPatches(); ++p )
2824// result.addPiecePointer(new gsShellStressFunction<d, T, bending>(m_patches,deformed,m_materialMatrices,p,type,m_assembler));
2825
2826// gsField<T> stressField(m_patches,result, true);
2827// return stressField;
2828
2829// }
2830
2831template<short_t d, class T, bool bending>
2832void gsThinShellAssembler<d, T, bending>::projectL2_into(const gsFunction<T> & /*fun*/, gsMatrix<T>& /*result*/)
2833{
2834 // /// todo: make a projection with BCs?
2835 // /// todo: test
2836 // // this->_getOptions();
2837
2838 // m_assembler.cleanUp();
2839 // GISMO_ENSURE(m_options.hasGroup("ExprAssembler"),"The option list does not contain options with the label 'ExprAssembler'!");
2840 // m_assembler.setOptions(m_options.getGroup("ExprAssembler"));
2841
2842 // geometryMap m_ori = m_assembler.getMap(m_patches);
2843
2844 // // Initialize stystem
2845 // m_assembler.initSystem();
2846
2847 // space m_space = m_assembler.trialSpace(0);
2848 // auto function = m_assembler.getCoeff(fun, m_ori);
2849 // // auto function = m_assembler.getCoeff(fun);
2850
2851 // // assemble system
2852 // m_assembler.assemble(m_space*m_space.tr()*meas(m_ori),m_space * function*meas(m_ori));
2853
2854 // gsSparseSolver<>::uPtr solver = gsSparseSolver<T>::get( m_options.askString("Solver","CGDiagonal") );
2855 // solver->compute(m_assembler.matrix());
2856 // result = solver->solve(m_assembler.rhs());
2857}
2858
2859template<short_t d, class T, bool bending>
2860void gsThinShellAssembler<d, T, bending>::projectL2_into(const gsFunction<T> & fun, gsMultiPatch<T>& mp)
2861{
2864 gsMatrix<T> tmp = projectL2(fun);
2865 mp = m_patches;
2866
2867 // Solution vector and solution variable
2868 space m_space = m_assembler.trialSpace(0);
2869 m_space.setup(m_bcs, dirichlet::l2Projection, m_continuity);
2870 const_cast<expr::gsFeSpace<T> & >(m_space).fixedPart() = m_ddofs;
2871
2872 solution m_solution = m_assembler.getSolution(m_space, tmp);
2873
2874 gsMatrix<T> cc;
2875 for ( size_t k =0; k!=mp.nPatches(); ++k) // Deform the geometry
2876 {
2877 // // extract deformed geometry
2878 m_solution.extract(cc, k);
2879 mp.patch(k).coefs() += cc; // defG points to mp_def, therefore updated
2880 }
2881}
2882
2883
2884template<short_t d, class T, bool bending>
2885gsMatrix<T> gsThinShellAssembler<d, T, bending>::projectL2(const gsFunction<T> & fun)
2886{
2889 gsMatrix<T> result;
2890 this->projectL2_into(fun,result);
2891 return result;
2892}
2893
2894template <short_t d, class T, bool bending>
2895T gsThinShellAssembler<d, T, bending>::deformationNorm(const gsMultiPatch<T> & deformed, const gsMultiPatch<T> & original)
2896{
2897 m_assembler.cleanUp();
2898 m_assembler.setOptions(m_options);
2899
2900 geometryMap ori = m_assembler.getMap(original);
2901 geometryMap def = m_assembler.getMap(deformed);
2902
2903 gsExprEvaluator<T> evaluator(m_assembler);
2904 T result = evaluator.integral(def.sqNorm() * meas(ori));
2905 return result;
2906}
2907
2908template<short_t d, class T, bool bending>
2909void gsThinShellAssembler<d, T, bending>::_ifcTest(const T tol)
2910{
2911 geometryMap m_ori = m_assembler.getMap(m_patches);
2912 gsExprEvaluator<T> ev(m_assembler);
2913
2914 m_inPlane.clear();
2915 m_outPlane.clear();
2916
2917 for (gsBoxTopology::const_iiterator it = m_patches.topology().iBegin(); it!=m_patches.topology().iEnd(); it++)
2918 {
2919 // G1 condition
2920 ev.integralInterface( (sn(m_ori.left()).normalized()-sn(m_ori.right()).normalized()).sqNorm() );
2921 ev.calcSqrt();
2922
2923 // // Continuous normal condition
2924 // ev.maxInterface( (sn(m_ori.left())-sn(m_ori.right())).norm() );
2925 // ev.calcSqrt();
2926
2927 if (ev.value() < tol)
2928 m_inPlane.push_back(*it);
2929 else
2930 m_outPlane.push_back(*it);
2931 }
2932}
2933
2934template<short_t d, class T, bool bending>
2935bool gsThinShellAssembler<d, T, bending>::_isInPlane(const boundaryInterface & /*ifc*/, const T tol)
2936{
2937 geometryMap m_ori = m_assembler.getMap(m_patches);
2938 gsExprEvaluator<T> ev(m_assembler);
2939
2940 // G1 condition
2941 ev.integralInterface( (sn(m_ori.left()).normalized()-sn(m_ori.right()).normalized()).sqNorm() );
2942 ev.calcSqrt();
2943
2944 // // Continuous normal condition
2945 // ev.maxInterface( (sn(m_ori.left())-sn(m_ori.right())).norm() );
2946 // ev.calcSqrt();
2947
2948 return (ev.value() < tol);
2949}
2950
2951
2952}// 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:2300
gismo::gsThinShellAssembler::assembleFoundationVector
ThinShellAssemblerStatus assembleFoundationVector(const gsFunctionSet< T > &deformed)
See gsThinShellAssemblerBase for details.
Definition gsThinShellAssembler.hpp:2270
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:2520
gismo::gsThinShellAssembler::constructDisplacement
gsMultiPatch< T > constructDisplacement(const gsMatrix< T > &solVector) const
See gsThinShellAssemblerBase for details.
Definition gsThinShellAssembler.hpp:2679
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:2895
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:2763
gismo::gsThinShellAssembler::assembleVector
ThinShellAssemblerStatus assembleVector(const gsFunctionSet< T > &deformed, const bool homogenize=true)
See gsThinShellAssemblerBase for details.
Definition gsThinShellAssembler.hpp:1980
gismo::gsThinShellAssembler::plotSolution
void plotSolution(std::string string, const gsMatrix< T > &solVector)
See gsThinShellAssemblerBase for details.
Definition gsThinShellAssembler.hpp:2563
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:2506
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:2193
gismo::gsThinShellAssembler::getElasticEnergy
T getElasticEnergy(const gsFunctionSet< T > &deformed)
See gsThinShellAssemblerBase for details.
Definition gsThinShellAssembler.hpp:2539
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:2795
gismo::gsThinShellAssembler::fullSolutionVector
gsMatrix< T > fullSolutionVector(const gsMatrix< T > &vector) const
See gsThinShellAssemblerBase for details.
Definition gsThinShellAssembler.hpp:2685
gismo::gsThinShellAssembler::assembleMatrix
ThinShellAssemblerStatus assembleMatrix(const gsFunctionSet< T > &deformed)
See gsThinShellAssemblerBase for details.
Definition gsThinShellAssembler.hpp:1703
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:2488
gismo::gsThinShellAssembler::projectL2_into
void projectL2_into(const gsFunction< T > &fun, gsMatrix< T > &result)
See gsThinShellAssemblerBase for details.
Definition gsThinShellAssembler.hpp:2832
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:2623
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:2697
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:2885
gismo::gsThinShellAssembler::_constructSolution
gsMultiPatch< T > _constructSolution(const gsMatrix< T > &solVector, const gsMultiPatch< T > &undeformed) const
See gsThinShellAssemblerBase for details.
Definition gsThinShellAssembler.hpp:2477
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:2131
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:2733
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