G+Smo  25.01.0
Geometry + Simulation Modules
 
Loading...
Searching...
No Matches
gsThinShellFunctions.hpp
Go to the documentation of this file.
1
16#pragma once
17
18#include <gsAssembler/gsExprEvaluator.h>
19
20namespace gismo
21{
22template <class T>
23void gsShellStressFunction<T>::eval_into(const gsMatrix<T> & u, gsMatrix<T> & result) const
24{
25 gsMatrix<T> Z(1,1);
26 Z.setZero();
27
28 result.setZero(targetDim(),u.cols());
29 gsMatrix<T> z(1,1);
30 z.setZero();
31
32 gsExprEvaluator<T> ev;
33
34 geometryMap m_ori = ev.getMap(*m_patches);
35 geometryMap m_def = ev.getMap(*m_defpatches);
36
37 // Initialize stystem
38 // m_assembler.initSystem(false);
39
40 gsMaterialMatrixIntegrate<T,MaterialOutput::VectorN> m_N0(m_materialMatrices,m_defpatches);
41 variable N0 = ev.getVariable(m_N0);
42 gsMaterialMatrixIntegrate<T,MaterialOutput::VectorM> m_N1(m_materialMatrices,m_defpatches);
43 variable N1 = ev.getVariable(m_N1);
44
45 gsMaterialMatrixIntegrate<T,MaterialOutput::CauchyVectorN> m_Nc0(m_materialMatrices,m_defpatches);
46 variable Nc0 = ev.getVariable(m_Nc0);
47 gsMaterialMatrixIntegrate<T,MaterialOutput::CauchyVectorM> m_Nc1(m_materialMatrices,m_defpatches);
48 variable Nc1 = ev.getVariable(m_Nc1);
49
50 gsMaterialMatrixEval<T,MaterialOutput::StressN> m_S0(m_materialMatrices,m_defpatches,z);
51 variable S0 = ev.getVariable(m_S0);
52 gsMaterialMatrixEval<T,MaterialOutput::StressM> m_S1(m_materialMatrices,m_defpatches,z);
53 variable S1 = ev.getVariable(m_S1);
54
55 gsMaterialMatrixEval<T,MaterialOutput::CauchyStressN> m_Sc0(m_materialMatrices,m_defpatches,z);
56 variable Sc0 = ev.getVariable(m_Sc0);
57 gsMaterialMatrixEval<T,MaterialOutput::CauchyStressM> m_Sc1(m_materialMatrices,m_defpatches,z);
58 variable Sc1 = ev.getVariable(m_Sc1);
59
60 gsMaterialMatrixEval<T,MaterialOutput::PStress> m_Sp(m_materialMatrices,m_defpatches,z);
61 variable Sp = ev.getVariable(m_Sp);
62 gsMaterialMatrixEval<T,MaterialOutput::PStressDir> m_pstressdir(m_materialMatrices,m_defpatches,z);
63 variable pstressdir = ev.getVariable(m_pstressdir);
64
65 gsMaterialMatrixEval<T,MaterialOutput::PCauchyStressN> m_Sp0(m_materialMatrices,m_defpatches,z);
66 variable Sp0 = ev.getVariable(m_Sp0);
67 gsMaterialMatrixEval<T,MaterialOutput::PCauchyStressM> m_Sp1(m_materialMatrices,m_defpatches,z);
68 variable Sp1 = ev.getVariable(m_Sp1);
69
70 gsMaterialMatrixEval<T,MaterialOutput::PStrainN> m_Ep0(m_materialMatrices,m_defpatches,z);
71 variable Ep0 = ev.getVariable(m_Ep0);
72 gsMaterialMatrixEval<T,MaterialOutput::PStrainM> m_Ep1(m_materialMatrices,m_defpatches,z);
73 variable Ep1 = ev.getVariable(m_Ep1);
74
75 gsMaterialMatrixEval<real_t,MaterialOutput::TensionField> m_TF(m_materialMatrices,m_defpatches,z);
76 variable TF = ev.getVariable(m_TF);
77
78 gsMaterialMatrixEval<T,MaterialOutput::Stretch> m_lambda(m_materialMatrices,m_defpatches,z);
79 variable lambda = ev.getVariable(m_lambda);
80 gsMaterialMatrixEval<T,MaterialOutput::StretchDir> m_lambdadir(m_materialMatrices,m_defpatches,z);
81 variable lambdadir = ev.getVariable(m_lambdadir);
82
83 gsFunctionExpr<> mult2t("1","0","0","0","1","0","0","0","2",2);
84 variable m_m2 = ev.getVariable(mult2t);
85 gsFunctionExpr<> mult12t("1","0","0","0","1","0","0","0","0.5",2);
86 variable m_m12 = ev.getVariable(mult12t);
87
88 auto That = cartcon(m_ori);
89 auto Ttilde = cartcov(m_ori);
90 auto Ttilde_def = cartcov(m_def);
91 // auto Tmat = cartcov(m_def);
92 auto E_m = 0.5 * ( flat(jac(m_def).tr()*jac(m_def)) - flat(jac(m_ori).tr()* jac(m_ori)) ) * reshape(m_m12,3,3) * That.tr();
93 auto E_f = ( expr::deriv2(m_ori,sn(m_ori).normalized().tr()) - expr::deriv2(m_def,sn(m_def).normalized().tr()) ) * reshape(m_m2,3,3) * That.tr();
94
95 auto N_m = N0.tr() * Ttilde.tr();
96 auto N_f = N1.tr() * Ttilde.tr();
97
98 auto Nc_m = Nc0.tr() * Ttilde_def.tr();
99 auto Nc_f = Nc1.tr() * Ttilde_def.tr();
100
101 auto S_m = S0.tr() * Ttilde.tr();
102 auto S_f = S1.tr() * Ttilde.tr();
103
104 auto Sc_m = Sc0.tr() * Ttilde_def.tr();
105 auto Sc_f = Sc1.tr() * Ttilde_def.tr();
106
107 gsMatrix<T> tmp;
108
109 switch (m_stress_type)
110 {
111 case stress_type::displacement :
112 for (index_t k = 0; k != u.cols(); ++k)
113 result.col(k) = ev.eval(m_def-m_ori,u.col(k),m_patchID);
114 break;
115
116 case stress_type::membrane :
117 for (index_t k = 0; k != u.cols(); ++k)
118 result.col(k) = (ev.eval(Sc_m,u.col(k),m_patchID)).transpose();
119 break;
120
121 case stress_type::flexural :
122 for (index_t k = 0; k != u.cols(); ++k)
123 result.col(k) = (ev.eval(Sc_f,u.col(k),m_patchID)).transpose();
124 break;
125
126 case stress_type::membrane_PK2 :
127 for (index_t k = 0; k != u.cols(); ++k)
128 result.col(k) = (ev.eval(S_m,u.col(k),m_patchID)).transpose();
129 break;
130
131 case stress_type::flexural_PK2 :
132 for (index_t k = 0; k != u.cols(); ++k)
133 result.col(k) = (ev.eval(S_f,u.col(k),m_patchID)).transpose();
134 break;
135
136 case stress_type::membrane_force :
137 for (index_t k = 0; k != u.cols(); ++k)
138 result.col(k) = (ev.eval(Nc_m,u.col(k),m_patchID)).transpose();
139 break;
140
141 case stress_type::flexural_moment :
142 for (index_t k = 0; k != u.cols(); ++k)
143 result.col(k) = (ev.eval(Nc_f,u.col(k),m_patchID)).transpose();
144 break;
145
146 case stress_type::membrane_force_PK2 :
147 for (index_t k = 0; k != u.cols(); ++k)
148 result.col(k) = (ev.eval(N_m,u.col(k),m_patchID)).transpose();
149 break;
150
151 case stress_type::flexural_moment_PK2 :
152 for (index_t k = 0; k != u.cols(); ++k)
153 result.col(k) = (ev.eval(N_f,u.col(k),m_patchID)).transpose();
154 break;
155
156 // -------------------------------------
157 case stress_type::von_mises :
158 for (index_t k = 0; k != u.cols(); ++k)
159 {
160 gsMatrix<> S;
161 gsMatrix<> Sm = (ev.eval(S_m,u.col(k),m_patchID)).transpose();
162 gsMatrix<> Sf = (ev.eval(S_f,u.col(k),m_patchID)).transpose();
163 S = Sm + Sf;
164 result(0,k) = math::sqrt(S(0,0)*S(0,0)+S(1,0)*S(1,0)-S(0,0)*S(1,0)+3*S(2,0)*S(2,0)); // ASSUMES PLANE STRESS
165 S = Sm - Sf;
166 result(1,k) = math::sqrt(S(0,0)*S(0,0)+S(1,0)*S(1,0)-S(0,0)*S(1,0)+3*S(2,0)*S(2,0)); // ASSUMES PLANE STRESS
167 }
168 break;
169
170 case stress_type::von_mises_membrane :
171 for (index_t k = 0; k != u.cols(); ++k)
172 {
173 gsMatrix<> S = (ev.eval(S_m,u.col(k),m_patchID)).transpose();
174 result(0,k) = math::sqrt(S(0,0)*S(0,0)+S(1,0)*S(1,0)-S(0,0)*S(1,0)+3*S(2,0)*S(2,0)); // ASSUMES PLANE STRESS
175 }
176 break;
177
178 case stress_type::von_mises_flexural :
179 for (index_t k = 0; k != u.cols(); ++k)
180 {
181 gsMatrix<> S = (ev.eval(S_f,u.col(k),m_patchID)).transpose();
182 result(0,k) = math::sqrt(S(0,0)*S(0,0)+S(1,0)*S(1,0)-S(0,0)*S(1,0)+3*S(2,0)*S(2,0)); // ASSUMES PLANE STRESS
183 }
184 break;
185 // -------------------------------------
186
187 case stress_type::membrane_strain :
188 for (index_t k = 0; k != u.cols(); ++k)
189 result.col(k) = ev.eval(E_m.tr(),u.col(k),m_patchID);
190 break;
191 case stress_type::flexural_strain :
192 for (index_t k = 0; k != u.cols(); ++k)
193 result.col(k) = ev.eval(E_f.tr(),u.col(k),m_patchID);
194 break;
195 case stress_type::principal_membrane_strain:
196 for (index_t k = 0; k != u.cols(); ++k)
197 result.col(k) = ev.eval(Ep0,u.col(k),m_patchID);
198 break;
199 case stress_type::principal_flexural_strain :
200 for (index_t k = 0; k != u.cols(); ++k)
201 result.col(k) = ev.eval(Ep1,u.col(k),m_patchID);
202 break;
203 case stress_type::principal_stretch :
204 for (index_t k = 0; k != u.cols(); ++k)
205 result.col(k) = ev.eval(lambda,u.col(k),m_patchID);
206 break;
207 case stress_type::principal_stress :
208 for (index_t k = 0; k != u.cols(); ++k)
209 result.col(k) = ev.eval(Sp,u.col(k));
210 break;
211 case stress_type::principal_stress_membrane :
212 for (index_t k = 0; k != u.cols(); ++k)
213 result.col(k) = ev.eval(Sp0,u.col(k),m_patchID);
214 break;
215 case stress_type::principal_stress_flexural :
216 for (index_t k = 0; k != u.cols(); ++k)
217 result.col(k) = ev.eval(Sp1,u.col(k),m_patchID);
218 break;
219 case stress_type::principal_stretch_dir1 :
220 for (index_t k = 0; k != u.cols(); ++k)
221 {
222 tmp = ev.eval(lambdadir,u.col(k),m_patchID);
223 result.col(k) = tmp.reshape(3,3).col(0);
224 }
225 break;
226 case stress_type::principal_stretch_dir2 :
227 for (index_t k = 0; k != u.cols(); ++k)
228 {
229 tmp = ev.eval(lambdadir,u.col(k),m_patchID);
230 result.col(k) = tmp.reshape(3,3).col(1);
231 }
232 break;
233 case stress_type::principal_stretch_dir3 :
234 for (index_t k = 0; k != u.cols(); ++k)
235 {
236 tmp = ev.eval(lambdadir,u.col(k),m_patchID);
237 result.col(k) = tmp.reshape(3,3).col(2);
238 }
239 break;
240 case stress_type::principal_stress_dir1 :
241 for (index_t k = 0; k != u.cols(); ++k)
242 {
243 tmp = ev.eval(pstressdir,u.col(k));
244 result.col(k) = tmp.reshape(3,3).col(0);
245 }
246 break;
247 case stress_type::principal_stress_dir2 :
248 for (index_t k = 0; k != u.cols(); ++k)
249 {
250 tmp = ev.eval(pstressdir,u.col(k));
251 result.col(k) = tmp.reshape(3,3).col(1);
252 }
253 break;
254 case stress_type::principal_stress_dir3 :
255 for (index_t k = 0; k != u.cols(); ++k)
256 {
257 tmp = ev.eval(pstressdir,u.col(k));
258 result.col(k) = tmp.reshape(3,3).col(2);
259 }
260 break;
261 case stress_type::tension_field :
262 for (index_t k = 0; k != u.cols(); ++k)
263 {
264 result.col(k) = ev.eval(TF,u.col(k),m_patchID);
265 }
266 break;
267 case stress_type::total :
268 gsWarn<<"Stress type 'total' not implemented\n";
269 break;
270 }
271}
272
273} // namespace gismo ends
gismo::gsExprEvaluator
Generic evaluator of isogeometric expressions.
Definition gsExprEvaluator.h:39
gismo::gsExprEvaluator::eval
void eval(const expr::_expr< E > &expr, gsGridIterator< T, mode, d > &git, const index_t patchInd=0)
gismo::gsExprEvaluator::getVariable
variable getVariable(const gsFunctionSet< T > &func, index_t dim=1)
Registers func as a variable and returns a handle to it.
Definition gsExprEvaluator.h:124
gismo::gsExprEvaluator::getMap
geometryMap getMap(const gsMultiPatch< T > &mp)
Registers mp as an isogeometric geometry map and return a handle to it.
Definition gsExprEvaluator.h:116
gismo::gsFunctionExpr
Class defining a multivariate (real or vector) function given by a string mathematical expression.
Definition gsFunctionExpr.h:52
gismo::gsMatrix
A matrix with arbitrary coefficient type and fixed or dynamic size.
Definition gsMatrix.h:41
gismo::gsMatrix::reshape
gsAsMatrix< T, Dynamic, Dynamic > reshape(index_t n, index_t m)
Returns the matrix resized to n x m matrix (data is not copied) This function assumes that the matrix...
Definition gsMatrix.h:221
gismo::gsShellStressFunction::eval_into
virtual void eval_into(const gsMatrix< T > &u, gsMatrix< T > &result) const
Each column of the input matrix (u) corresponds to one evaluation point. Each column of the output ma...
Definition gsThinShellFunctions.hpp:23
index_t
#define index_t
Definition gsConfig.h:32
gsWarn
#define gsWarn
Definition gsDebug.h:50
gsExprEvaluator.h
Generic expressions evaluator.
gismo
The G+Smo namespace, containing all definitions for the library.
gismo::stress_type::principal_stress_dir1
@ principal_stress_dir1
principal stretch directions
Definition gsThinShellFunctions.h:64
gismo::stress_type::membrane_force
@ membrane_force
compute membrane Cauchy stresses
Definition gsThinShellFunctions.h:45
gismo::stress_type::von_mises_flexural
@ von_mises_flexural
compute only von Mises stress - membrane stresses
Definition gsThinShellFunctions.h:43
gismo::stress_type::membrane
@ membrane
compute only von Mises stress - flexural stresses
Definition gsThinShellFunctions.h:44
gismo::stress_type::membrane_force_PK2
@ membrane_force_PK2
compute membrane PK2 stresses
Definition gsThinShellFunctions.h:47
gismo::stress_type::flexural_moment_PK2
@ flexural_moment_PK2
compute flexural Cauchy stresses integrated over the thickness
Definition gsThinShellFunctions.h:51
gismo::stress_type::von_mises_membrane
@ von_mises_membrane
compute only von Mises stress
Definition gsThinShellFunctions.h:42
gismo::stress_type::membrane_PK2
@ membrane_PK2
compute membrane Cauchy stresses integrated over the thickness
Definition gsThinShellFunctions.h:46
gismo::stress_type::principal_stress
@ principal_stress
principal stretches
Definition gsThinShellFunctions.h:58
gismo::stress_type::tension_field
@ tension_field
principal stress directions
Definition gsThinShellFunctions.h:67
gismo::stress_type::principal_stress_membrane
@ principal_stress_membrane
principal stress membrane
Definition gsThinShellFunctions.h:59
gismo::stress_type::principal_stress_dir3
@ principal_stress_dir3
principal stress directions
Definition gsThinShellFunctions.h:66
gismo::stress_type::principal_stretch_dir2
@ principal_stretch_dir2
principal stretch directions
Definition gsThinShellFunctions.h:62
gismo::stress_type::principal_stress_flexural
@ principal_stress_flexural
principal stress membrane
Definition gsThinShellFunctions.h:60
gismo::stress_type::principal_stretch_dir1
@ principal_stretch_dir1
principal stress bending
Definition gsThinShellFunctions.h:61
gismo::stress_type::flexural_moment
@ flexural_moment
compute flexural PK2 stresses
Definition gsThinShellFunctions.h:50
gismo::stress_type::flexural
@ flexural
compute membrane PK2 stresses integrated over the thickness
Definition gsThinShellFunctions.h:48
gismo::stress_type::principal_stress_dir2
@ principal_stress_dir2
principal stress directions
Definition gsThinShellFunctions.h:65
gismo::stress_type::principal_stretch_dir3
@ principal_stretch_dir3
principal stretch directions
Definition gsThinShellFunctions.h:63
gismo::stress_type::total
@ total
compute flexural PK2 stresses integrated over the thickness
Definition gsThinShellFunctions.h:52
gismo::stress_type::flexural_PK2
@ flexural_PK2
compute flexural Cauchy stresses
Definition gsThinShellFunctions.h:49