G+Smo  25.01.0
Geometry + Simulation Modules
 
Loading...
Searching...
No Matches
gsMaterialMatrixTFT.hpp
Go to the documentation of this file.
1
17#pragma once
18
19#include <gsKLShell/src/gsMaterialMatrixUtils.h>
20#include <gsCore/gsFunction.h>
21
22template <typename T> T mod(const T x, const T y);
23template <typename T> class objective;
24
25namespace gismo
26{
27
28template <short_t dim, class T, bool linear >
29std::ostream & gsMaterialMatrixTFT<dim,T,linear>::print(std::ostream &os) const
30{
31 os <<"---------------------------------------------------------------------\n"
32 <<"---------------------Tension-Field Theory Material-------------------\n"
33 <<"---------------------------------------------------------------------\n\n";
34
35 m_materialMat->print(os);
36
37 os <<"---------------------------------------------------------------------\n\n";
38 return os;
39}
40
41template <short_t dim, class T, bool linear >
42gsMatrix<T> gsMaterialMatrixTFT<dim,T,linear>::eval3D_matrix(const index_t patch, const gsMatrix<T> & u, const gsMatrix<T> & z, enum MaterialOutput /*out*/) const
43{
44 // GISMO_ASSERT(out==MaterialOutput::MatrixA,"Tension Field Theory only works for membrane models, hence only outputs the A matrix");
45 // Input: u in-plane points
46 // z matrix with, per point, a column with z integration points
47 // Output: (n=u.cols(), m=z.rows())
48 // [(u1,z1) (u2,z1) .. (un,z1), (u1,z2) .. (un,z2), .., (u1,zm) .. (un,zm)]
49
50 // return this->_eval3D_matrix_impl<true>(patch,u,z);
51 return this->_eval3D_matrix_impl<linear>(patch,u,z);
52}
53
54template <short_t dim, class T, bool linear >
55gsMatrix<T> gsMaterialMatrixTFT<dim,T,linear>::eval3D_vector(const index_t patch, const gsMatrix<T> & u, const gsMatrix<T> & z, enum MaterialOutput out) const
56{
57 return this->eval3D_stress(patch,u,z,out);
58}
59
60template <short_t dim, class T, bool linear >
61gsMatrix<T> gsMaterialMatrixTFT<dim,T,linear>::eval3D_pstress(const index_t patch, const gsMatrix<T> & u, const gsMatrix<T>& z, enum MaterialOutput out) const
62{
63 gsMatrix<T> TF = this->_compute_TF(patch,u,z);
64 gsMatrix<T> result = m_materialMat->eval3D_pstress(patch,u,z,out);
65 index_t colIdx;
66 for (index_t k=0; k!=u.cols(); k++)
67 {
68 for( index_t j=0; j < z.rows(); ++j ) // through-thickness points
69 {
70 colIdx = j*u.cols()+k;
71 if (TF(0,colIdx) == 1 || TF(0,colIdx) == 0) // taut & wrinkled
72 {
73 // do nothing
74 }
75 else if (TF(0,colIdx) == -1) // slack
76 {
77 // Set to zero
78 // result.col(colIdx).setZero();
79 result.col(colIdx) *= m_options.getReal("SlackMultiplier");
80 }
81 else
82 GISMO_ERROR("Tension field data not understood: "<<TF(0,colIdx));
83 }
84 }
85 return result;
86}
87
88template <short_t dim, class T, bool linear >
89gsMatrix<T> gsMaterialMatrixTFT<dim,T,linear>::eval3D_pstrain(const index_t patch, const gsMatrix<T> & u, const gsMatrix<T>& z) const
90{
91 gsMatrix<T> TF = this->_compute_TF(patch,u,z);
92 gsMatrix<T> result = m_materialMat->eval3D_pstrain(patch,u,z);
93 index_t colIdx;
94 for (index_t k=0; k!=u.cols(); k++)
95 {
96 for( index_t j=0; j < z.rows(); ++j ) // through-thickness points
97 {
98 colIdx = j*u.cols()+k;
99 if (TF(0,colIdx) == 1 || TF(0,colIdx) == 0) // taut & wrinkled
100 {
101 // do nothing
102 }
103 else if (TF(0,colIdx) == -1) // slack
104 {
105 // Set to zero
106 // result.col(colIdx).setZero();
107 result.col(colIdx) *= m_options.getReal("SlackMultiplier");
108 }
109 else
110 GISMO_ERROR("Tension field data not understood: "<<TF(0,colIdx));
111 }
112 }
113 return result;
114}
115
116template <short_t dim, class T, bool linear >
117gsMatrix<T> gsMaterialMatrixTFT<dim,T,linear>::eval3D_strain(const index_t patch, const gsMatrix<T> & u, const gsMatrix<T> & z) const
118{
119 gsMatrix<T> result = m_materialMat->eval3D_strain(patch,u,z);
120 gsMatrix<T> TF = this->_compute_TF(patch,u,z);
121 index_t colIdx;
122 T theta;
123 for (index_t k=0; k!=u.cols(); k++)
124 {
125 for( index_t j=0; j < z.rows(); ++j ) // through-thickness points
126 {
127 colIdx = j*u.cols()+k;
128 if (TF(0,colIdx) == 1) // taut
129 {
130 // do nothing
131 }
132 else if (TF(0,colIdx) == -1) // slack
133 {
134 // Set to zero
135 // result.col(colIdx).setZero();
136 result.col(colIdx) *= m_options.getReal("SlackMultiplier");
137 }
138 else if (TF(0,colIdx) == 0) // wrinkled
139 {
140 gsMatrix<T> C = m_materialMat->eval3D_matrix(patch,u.col(k),z(j,k),MaterialOutput::MatrixA);
141 gsMatrix<T> N = m_materialMat->eval3D_stress(patch,u.col(k),z(j,k),MaterialOutput::VectorN);
142 gsMatrix<T> E = m_materialMat->eval3D_strain(patch,u.col(k),z(j,k));
143 gsMatrix<T> thetas = eval_theta(C,N,E);
144 theta = thetas(0,0);
145 result.col(colIdx) = this->_compute_E(theta,C.reshape(3,3),N,E);
146 }
147 else
148 GISMO_ERROR("Tension field data not understood: "<<TF(0,colIdx));
149 }
150 }
151 return result;
152}
153
154template <short_t dim, class T, bool linear >
155gsMatrix<T> gsMaterialMatrixTFT<dim,T,linear>::eval3D_stress(const index_t patch, const gsMatrix<T> & u, const gsMatrix<T> & z, enum MaterialOutput /*out*/) const
156{
157 // Input: u in-plane points
158 // z matrix with, per point, a column with z integration points
159 // Output: (n=u.cols(), m=z.rows())
160 // [(u1,z1) (u2,z1) .. (un,z1), (u1,z2) .. (un,z2), .., (u1,zm) .. (un,zm)]
161
162 gsMatrix<T> result = m_materialMat->eval3D_stress(patch,u,z,MaterialOutput::VectorN);
163 gsMatrix<T> TF = this->_compute_TF(patch,u,z);
164 index_t colIdx;
165 T theta;
166 for (index_t k=0; k!=u.cols(); k++)
167 {
168 for( index_t j=0; j < z.rows(); ++j ) // through-thickness points
169 {
170 colIdx = j*u.cols()+k;
171 if (TF(0,colIdx) == 1) // taut
172 {
173 // do nothing
174 }
175 else if (TF(0,colIdx) == -1) // slack
176 {
177 // Set to zero
178 // result.col(colIdx).setZero();
179 result.col(colIdx) *= m_options.getReal("SlackMultiplier");
180 }
181 else if (TF(0,colIdx) == 0) // wrinkled
182 {
183 gsMatrix<T> C = m_materialMat->eval3D_matrix(patch,u.col(k),z(j,k),MaterialOutput::MatrixA);
184 gsAsMatrix<T,Dynamic,Dynamic> N = result.reshapeCol(colIdx,3,1);
185 gsMatrix<T> E = m_materialMat->eval3D_strain(patch,u.col(k),z(j,k));
186
187 gsMatrix<T> thetas = eval_theta(C,N,E);
188 theta = thetas(0,0);
189 result.col(colIdx) = this->_compute_S(theta,C.reshape(3,3),N);
190 }
191 else
192 GISMO_ERROR("Tension field data not understood: "<<TF(0,colIdx));
193 }
194 }
195 return result;
196}
197
198template <short_t dim, class T, bool linear >
199gsMatrix<T> gsMaterialMatrixTFT<dim,T,linear>::eval3D_CauchyStress(const index_t patch, const gsMatrix<T> & u, const gsMatrix<T> & z, enum MaterialOutput /*out*/) const
200{
201 // Input: u in-plane points
202 // z matrix with, per point, a column with z integration points
203 // Output: (n=u.cols(), m=z.rows())
204 // [(u1,z1) (u2,z1) .. (un,z1), (u1,z2) .. (un,z2), .., (u1,zm) .. (un,zm)]
205
206 gsMatrix<T> result = m_materialMat->eval3D_CauchyStress(patch,u,z,MaterialOutput::CauchyStressN);
207 gsMatrix<T> TF = this->_compute_TF(patch,u,z);
208 index_t colIdx;
209 T theta;
210
211 this->_computePoints(patch,u);
212 for (index_t k=0; k!=u.cols(); k++)
213 {
214 for( index_t j=0; j < z.rows(); ++j ) // through-thickness points
215 {
216 colIdx = j*u.cols()+k;
217 if (TF(0,colIdx) == 1) // taut
218 {
219 // do nothing
220 }
221 else if (TF(0,colIdx) == -1) // slack
222 {
223 // Set to zero
224 // result.col(colIdx).setZero();
225 result.col(colIdx) *= m_options.getReal("SlackMultiplier");
226 }
227 else if (TF(0,colIdx) == 0) // wrinkled
228 {
229 this->_getMetric(k, z(j, k) * m_data.mine().m_Tmat(0, k)); // on point i, on height z(0,j)
230 gsMatrix<T> C = m_materialMat->eval3D_matrix(patch,u.col(k),z(j,k),MaterialOutput::MatrixA);
231 gsMatrix<T> N = m_materialMat->eval3D_stress(patch,u.col(k),z(j,k),MaterialOutput::VectorN);;
232 gsMatrix<T> E = m_materialMat->eval3D_strain(patch,u.col(k),z(j,k));
233 gsMatrix<T> thetas = eval_theta(C,N,E);
234 theta = thetas(0,0);
235 gsMatrix<T> S = this->_compute_S(theta,C.reshape(3,3),N);
236 T detF = math::sqrt(m_data.mine().m_J0_sq)*1.0;
237 result.col(colIdx) = S / math::sqrt(detF);
238 }
239 else
240 GISMO_ERROR("Tension field data not understood: "<<TF(0,colIdx));
241 }
242 }
243 return result;
244}
245
246// template <short_t dim, class T, bool linear >
247// gsMatrix<T> gsMaterialMatrixTFT<dim,T,linear>::eval3D_CauchyPStress(const index_t patch, const gsMatrix<T> & u, const gsMatrix<T> & z, enum MaterialOutput out) const
248// {
249// // Input: u in-plane points
250// // z matrix with, per point, a column with z integration points
251// // Output: (n=u.cols(), m=z.rows())
252// // [(u1,z1) (u2,z1) .. (un,z1), (u1,z2) .. (un,z2), .., (u1,zm) .. (un,zm)]
253
254// gsMatrix<T> result = m_materialMat->eval3D_CauchyPStress(patch,u,z,MaterialOutput::PCauchyStressN);
255// gsMatrix<T> TF = this->_compute_TF(patch,u,z);
256// index_t colIdx;
257// T theta;
258// std::pair<gsVector<T>,gsMatrix<T>> res;
259
260// this->_computePoints(patch,u);
261// for (index_t k=0; k!=u.cols(); k++)
262// {
263// for( index_t j=0; j < z.rows(); ++j ) // through-thickness points
264// {
265// colIdx = j*u.cols()+k;
266// if (TF(0,colIdx) == 1) // taut
267// {
268// // do nothing
269// }
270// else if (TF(0,colIdx) == -1) // slack
271// {
272// // Set to zero
273// // result.col(colIdx).setZero();
274// result.col(colIdx) *= m_options.getReal("SlackMultiplier");
275// }
276// else if (TF(0,colIdx) == 0) // wrinkled
277// {
278// this->_getMetric(k, z(j, k) * m_data.mine().m_Tmat(0, k)); // on point i, on height z(0,j)
279// gsMatrix<T> C = m_materialMat->eval3D_matrix(patch,u.col(k),z(j,k),MaterialOutput::MatrixA);
280// gsMatrix<T> N = m_materialMat->eval3D_stress(patch,u.col(k),z(j,k),MaterialOutput::VectorN);;
281// gsMatrix<T> E = m_materialMat->eval3D_strain(patch,u.col(k),z(j,k),MaterialOutput::StrainN);
282// gsMatrix<T> thetas = eval_theta(C,N,E);
283// theta = thetas(0,0);
284// gsMatrix<T> S = this->_compute_S(theta,C.reshape(3,3),N);
285// res = this->_evalPStress(S);
286
287// T detF = math::sqrt(m_data.mine().m_J0_sq)*1.0;
288// result.col(colIdx) = res.first / math::sqrt(detF);
289// }
290// else
291// GISMO_ERROR("Tension field data not understood: "<<TF(0,colIdx));
292// }
293// }
294// return result;
295// }
296
297template <short_t dim, class T, bool linear >
298gsMatrix<T> gsMaterialMatrixTFT<dim,T,linear>::eval3D_tensionfield(const index_t patch, const gsMatrix<T> & u, const gsMatrix<T>& z, enum MaterialOutput /*out*/) const
299{
300 return this->_compute_TF(patch,u,z);
301}
302
303template <short_t dim, class T, bool linear >
304gsMatrix<T> gsMaterialMatrixTFT<dim,T,linear>::eval3D_theta(const index_t patch, const gsMatrix<T> & u, const gsMatrix<T>& z, enum MaterialOutput /*out*/) const
305{
306 gsMatrix<T> result(1,u.cols());
307 result.setZero();
308 gsMatrix<T> TF = this->_compute_TF(patch,u,z);
309 index_t colIdx;
310 for (index_t k=0; k!=u.cols(); k++)
311 {
312 for( index_t j=0; j < z.rows(); ++j ) // through-thickness points
313 {
314 colIdx = j*u.cols()+k;
315 if (TF(0,colIdx) == 1 || TF(0,colIdx) == -1) // taut
316 {
317 result(0,colIdx) = NAN;
318 }
319 else if (TF(0,colIdx) == 0) // wrinkled
320 {
321 gsMatrix<T> C = m_materialMat->eval3D_matrix(patch,u.col(k),z(j,k),MaterialOutput::MatrixA);
322 gsMatrix<T> N = m_materialMat->eval3D_stress(patch,u.col(k),z(j,k),MaterialOutput::VectorN);;
323 gsMatrix<T> E = m_materialMat->eval3D_strain(patch,u.col(k),z(j,k));
324
325 gsMatrix<T> thetas = eval_theta(C,N,E);
326 result(0,colIdx) = thetas(0,0);
327 }
328 else
329 GISMO_ERROR("Tension field data not understood: "<<TF(0,colIdx));
330 }
331 }
332 return result;
333}
334
335template <short_t dim, class T, bool linear >
336gsMatrix<T> gsMaterialMatrixTFT<dim,T,linear>::eval3D_gamma(const index_t patch, const gsMatrix<T> & u, const gsMatrix<T>& z, enum MaterialOutput /*out*/) const
337{
338 gsMatrix<T> result(1,u.cols());
339 result.setZero();
340 gsMatrix<T> TF = this->_compute_TF(patch,u,z);
341 index_t colIdx;
342 for (index_t k=0; k!=u.cols(); k++)
343 {
344 for( index_t j=0; j < z.rows(); ++j ) // through-thickness points
345 {
346 colIdx = j*u.cols()+k;
347 if (TF(0,colIdx) == 1 || TF(0,colIdx) == -1) // taut
348 {
349 result(0,colIdx) = NAN;
350 }
351 else if (TF(0,colIdx) == 0) // wrinkled
352 {
353 gsMatrix<T> C = m_materialMat->eval3D_matrix(patch,u.col(k),z(j,k),MaterialOutput::MatrixA);
354 gsMatrix<T> N = m_materialMat->eval3D_stress(patch,u.col(k),z(j,k),MaterialOutput::VectorN);;
355 gsMatrix<T> E = m_materialMat->eval3D_strain(patch,u.col(k),z(j,k));
356
357 gsMatrix<T> thetas = eval_theta(C,N,E);
358 result(0,colIdx) = _compute_gamma(thetas(0,0),C.reshape(3,3),N.reshape(3,1),E.reshape(3,1));
359 }
360 else
361 GISMO_ERROR("Tension field data not understood: "<<TF(0,colIdx));
362 }
363 }
364 return result;
365}
366
368template <short_t dim, class T, bool linear >
369template <bool _linear>
370typename std::enable_if<_linear, gsMatrix<T> >::type
371gsMaterialMatrixTFT<dim,T,linear>::_eval3D_matrix_impl(const index_t patch, const gsMatrix<T> & u, const gsMatrix<T>& z, enum MaterialOutput /*out*/) const
372{
373 gsMatrix<T> result = m_materialMat->eval3D_matrix(patch,u,z,MaterialOutput::MatrixA);
374 gsMatrix<T> TF = this->_compute_TF(patch,u,z);
375
376 index_t colIdx;
377 T theta;
378 for (index_t k=0; k!=u.cols(); k++)
379 {
380 for( index_t j=0; j < z.rows(); ++j ) // through-thickness points
381 {
382 colIdx = j*u.cols()+k;
383 if (TF(0,colIdx) == 1) // taut
384 {
385 // do nothing
386 }
387 else if (TF(0,colIdx) == -1) // slack
388 {
389 // Set to zero
390 // result.col(colIdx).setZero();
391 result.col(colIdx) *= m_options.getReal("SlackMultiplier");
392 }
393 else if (TF(0,colIdx) == 0) // wrinkled
394 {
395 gsMatrix<T> C = result.col(colIdx);
396 gsMatrix<T> N = m_materialMat->eval3D_stress(patch,u.col(k),z(j,k),MaterialOutput::VectorN);
397 gsMatrix<T> E = m_materialMat->eval3D_strain(patch,u.col(k),z(j,k));
398 gsMatrix<T> thetas = eval_theta(C,N,E);
399
400 theta = thetas(0,0);
401 result.col(colIdx) = this->_compute_C(theta,C.reshape(3,3),N.reshape(3,1)).reshape(9,1);
402 }
403 else
404 GISMO_ERROR("Tension field data not understood: "<<TF(0,colIdx));
405 }
406 }
407 return result;
408}
409
410template <short_t dim, class T, bool linear >
411template <bool _linear>
412typename std::enable_if<!_linear, gsMatrix<T> >::type
413gsMaterialMatrixTFT<dim,T,linear>::_eval3D_matrix_impl(const index_t patch, const gsMatrix<T> & u, const gsMatrix<T>& z, enum MaterialOutput out) const
414{
415 gsMatrix<T> result = m_materialMat->eval3D_matrix(patch,u,z,out);
416 gsMatrix<T> TF = this->_compute_TF(patch,u,z);
417
418 index_t colIdx;
419 T theta;
420 for (index_t k=0; k!=u.cols(); k++)
421 {
422 for( index_t j=0; j < z.rows(); ++j ) // through-thickness points
423 {
424 colIdx = j*u.cols()+k;
425 if (TF(0,colIdx) == 1) // taut
426 {
427 // do nothing
428 }
429 else if (TF(0,colIdx) == -1) // slack
430 {
431 // Set to zero
432 // result.col(colIdx).setZero();
433 result.col(colIdx) *= m_options.getReal("SlackMultiplier");
434 }
435 else if (TF(0,colIdx) == 0) // wrinkled
436 {
437 gsMatrix<T> C = result.col(colIdx);
438 gsMatrix<T> N = m_materialMat->eval3D_stress(patch,u.col(k),z(j,k),MaterialOutput::VectorN);
439 gsMatrix<T> E = m_materialMat->eval3D_strain(patch,u.col(k),z(j,k));
440 gsMatrix<T> thetas = eval_theta(C,N,E);
441
442 gsMatrix<T> dC = m_materialMat->eval3D_dmatrix(patch,u.col(k),z(j,k),MaterialOutput::MatrixA);
443 dC *= 2; // since eval3D_dmatrix gives dC/dC and dC/dE = 2*dC/dC
444 theta = thetas(0,0);
445 result.col(colIdx) = this->_compute_C(theta,C.reshape(3,3),N.reshape(3,1),dC.reshape(9,3)).reshape(9,1);
446 }
447 else
448 GISMO_ERROR("Tension field data not understood: "<<TF(0,colIdx));
449 }
450 }
451 return result;
452}
453
454template <short_t dim, class T, bool linear >
455gsMatrix<T> gsMaterialMatrixTFT<dim,T,linear>::eval_theta(const gsMatrix<T> & Cs, const gsMatrix<T> & Ns, const gsMatrix<T> & Es) const
456{
457 GISMO_ASSERT(Cs.cols()==Ns.cols(),"Number of C matrices and N vectors is different");
458 GISMO_ASSERT(Cs.cols()==Es.cols(),"Number of C matrices and E vectors is different");
459 gsMatrix<T> result(1,Cs.cols());
460 gsVector<T,3> n1_vec, n2_vec, n4_vec;
461
462 T theta = 0;
463
464 for (index_t k = 0; k!=Cs.cols(); k++)
465 {
466 gsMatrix<T> C = Cs.col(k);
467 C.resize(3,3);
468 gsMatrix<T> N = Ns.col(k);
469 N.resize(3,1);
470 gsMatrix<T> E = Es.col(k);
471 E.resize(3,1);
472
473 gsVector<T> theta_interval = _theta_interval(C, N, E);
474 objective<T> obj(C,N);
475
476 gsVector<T> zeros(1); zeros<<0;
477 gsVector<T> arg(1); //arg<<theta;
478 T f;
479 bool converged = false;
480 if (theta_interval(0)!=0 && theta_interval(1)!=0)
481 {
482 arg<<theta_interval(0);
483 converged = obj.findRootBrent(f,theta_interval(0),theta_interval(1),theta,1e-18);
484 if (!converged)
485 {
486 converged = obj.findRootBisection(f,theta_interval(0),theta_interval(1),theta,1e-6);
487 if (!converged)
488 {
489 obj.newtonRaphson(zeros,arg,false,1e-12,100,1);
490 theta = arg(0);
491 f = obj.eval(theta);
492 }
493 }
494 }
495
496 if (!_check_theta_full(theta,C,N,E) || (theta_interval(0)==0 && theta_interval(1)==0))
497 {
498 // Take the whole interval and evaluate 50 points
499 gsVector<T> theta_vec = gsVector<T>::LinSpaced(51,0,EIGEN_PI);
500 gsVector<T> f_vec(theta_vec.size());
501 // Compute f for each point
502 T t;
503 for (index_t i=0; i!=theta_vec.size(); i++)
504 {
505 t = theta_vec(i);
506 f_vec(i) = obj.eval(t);
507 }
508 // For each interval, store the intervals where the sign of f changes
509 std::vector<std::pair<T,T>> intervals; intervals.reserve(theta_vec.size());
510 for (index_t i=0; i!=f_vec.size()-1; i++)
511 if (f_vec(i)*f_vec(i+1) < 0)
512 intervals.push_back(std::make_pair(theta_vec(i),theta_vec(i+1)));
513
514 // Compute all intervals
515 std::vector<T> thetas(intervals.size());
516 for (size_t i = 0; i!=intervals.size(); i++)
517 {
518 arg<<intervals[i].first;
519 converged = obj.findRootBrent(f,intervals[i].first,intervals[i].second,theta,1e-18);
520 if (!converged)
521 {
522 converged = obj.findRootBisection(f,intervals[i].first,intervals[i].second,theta,1e-6);
523 if (!converged)
524 {
525 obj.newtonRaphson(zeros,arg,false,1e-12,100,1);
526 theta = arg(0);
527 f = obj.eval(theta);
528 }
529 }
530 thetas[i] = theta;
531 }
532
533 if ((f_vec(0)*f_vec(f_vec.size()-1) <= 0))
534 thetas.push_back(0.0);
535
536 std::vector<bool> full_check(thetas.size());
537 for (size_t i = 0; i!=thetas.size(); i++)
538 full_check[i] = _check_theta_full(thetas[i],C,N,E);
539
540 index_t full_check_sum = std::accumulate(full_check.begin(),full_check.end(),0);
541 if (full_check_sum!=0)
542 {
543 std::vector<bool>::iterator res = std::find_if(full_check.begin(), full_check.end(), [](bool i) { return i;});
544 index_t i = std::distance(full_check.begin(),res);
545 theta = thetas[i];
546 }
547 else
548 {
549 gsWarn<<"Everything failed??\n";
550 // gsDebugVar(E);
551 // gsDebugVar(N);
552 // gsDebugVar(C);
553 // gsVector<T> x = gsVector<T>::LinSpaced(1000,-EIGEN_PI,EIGEN_PI);
554 // for (index_t XX=0; XX!=x.size(); XX++)
555 // {
556 // T t = x(XX);
557 // f = obj.eval(t);
558
559 // n1 = math::cos(t);
560 // n2 = math::sin(t);
561 // m1 = -math::sin(t);
562 // m2 = math::cos(t);
563
564 // n1_vec<<n1*n1, n2*n2, 2*n1*n2;
565 // n2_vec<<m1*n1, m2*n2, m1*n2+m2*n1;
566 // gsVector<T,3> n4_vec; n4_vec<<m1*m1, m2*m2, 2*m1*m2;
567
568 // gamma = - ( n1_vec.transpose() * N ).value() / ( n1_vec.transpose() * C * n1_vec ).value();
569
570 // gsInfo<<x(XX)<<","<<f<<","<<gamma<<","<<(n1_vec.transpose() * N).value()<<","<<(n4_vec.transpose() * Es.col(k)).value()<<"\n";
571 // }
572 theta = 0;
573 }
574 }
575 result(0,k) = theta;
576 }
577 return result;
578}
579
580template <short_t dim, class T, bool linear >
581T gsMaterialMatrixTFT<dim,T,linear>::_compute_gamma(const T & theta, const gsMatrix<T> & C, const gsVector<T> & N, const gsVector<T> & /*E*/) const
582{
583 T n1 = math::cos(theta);
584 T n2 = math::sin(theta);
585 T m1 = -math::sin(theta);
586 T m2 = math::cos(theta);
587 gsVector<T,3> n1_vec; n1_vec<<n1*n1, n2*n2, 2*n1*n2;
588 gsVector<T,3> n2_vec; n2_vec<<m1*n1, m2*n2, m1*n2+m2*n1;
589 gsVector<T,3> n4_vec; n4_vec<<m1*m1, m2*m2, 2*m1*m2;
590
591 return - ( n1_vec.transpose() * N ).value() / ( n1_vec.transpose() * C * n1_vec ).value();
592}
593
594template <short_t dim, class T, bool linear >
595bool gsMaterialMatrixTFT<dim,T,linear>::_check_theta_full(const T & theta, const gsMatrix<T> & /*C*/, const gsVector<T> & N, const gsVector<T> & E) const
596{
597 T n1 = math::cos(theta);
598 T n2 = math::sin(theta);
599 T m1 = -math::sin(theta);
600 T m2 = math::cos(theta);
601 gsVector<T,3> n1_vec; n1_vec<<n1*n1, n2*n2, 2*n1*n2;
602 gsVector<T,3> n4_vec; n4_vec<<m1*m1, m2*m2, 2*m1*m2;
603
604 bool check = true;
605 check &= ((n1_vec.transpose() * N).value() < 0); // Li et al eq. 63 / eq. 55
606 check &= ((n4_vec.transpose() * E).value() > 0); // Li et al eq. 58
607 return check;
608}
609
610template <short_t dim, class T, bool linear >
611bool gsMaterialMatrixTFT<dim,T,linear>::_check_theta_gamma(const T & theta, const gsMatrix<T> & /*C*/, const gsVector<T> & N, const gsVector<T> & /*E*/) const
612{
613 T n1 = math::cos(theta);
614 T n2 = math::sin(theta);
615 gsVector<T,3> n1_vec; n1_vec<<n1*n1, n2*n2, 2*n1*n2;
616 return ((n1_vec.transpose() * N).value() < 0); // Li et al eq. 63 / eq. 55
617}
618
619template <short_t dim, class T, bool linear >
620gsVector<T> gsMaterialMatrixTFT<dim,T,linear>::_theta_interval(const gsMatrix<T> & C, const gsVector<T> & N, const gsVector<T> & E) const
621{
622 GISMO_ASSERT(C.rows()==C.cols(),"C must be a 3x3 square matrix");
623 GISMO_ASSERT(C.rows()==3,"C must be a 3x3 square matrix");
624 GISMO_ASSERT(N.rows()==3,"N must be a 3x1 vector");
625 GISMO_ASSERT(E.rows()==3,"E must be a 3x1 vector");
626 gsVector<T> result(2);
627
628 objective<T> obj(C,N);
629
630 // See Lu et al., Finite Element Analysis of Membrane Wrinkling, 2001, International Journal for numerical methods in engineering
631 T E11 = E(0);
632 T E22 = E(1);
633 T E12 = 0.5 * E(2);
634 T R_E = math::sqrt( math::pow((E11-E22)/2.,2) + E12*E12 );
635 if (R_E==0) {R_E = 1; gsDebugVar("??");};
636 T sin_E0 = E12/R_E;
637 T cos_E0 = (E22-E11)/(2*R_E);
638 std::complex<T> sin_Esqrt = E12*E12-E11*E22;
639 std::complex<T> sin_E1_C = -math::sqrt(sin_Esqrt)/R_E;
640 T sin_E1 = sin_E1_C.real();
641 T cos_E1 = -(E11+E22)/(2*R_E);
642 std::complex<T> sin_E2_C = math::sqrt(sin_Esqrt)/R_E;
643 T sin_E2 = sin_E2_C.real();
644 T cos_E2 = -(E11+E22)/(2*R_E);
645
646 T theta_0E = math::atan2(sin_E0,cos_E0);
647 T theta_1E = math::atan2(sin_E1,cos_E1);
648 T theta_2E = math::atan2(sin_E2,cos_E2);
649
650 T pi = 4*math::atan(1);
651 T N11 = N(0);
652 T N22 = N(1);
653 T N12 = N(2);
654 T R_N = math::sqrt( math::pow((N11-N22)/2.,2) + N12*N12 );
655 if (R_N==0) R_N = 1;
656 T sin_N0 = -N12/R_N;
657 T cos_N0 = (N11-N22)/(R_N);
658 std::complex<T> sin_Nsqrt = N12*N12-N11*N22;
659 // if (sin_Nsqrt < 0) gsDebugVar("Oops");
660 std::complex<T> sin_N1_C = math::sqrt(sin_Nsqrt)/R_N;
661 T sin_N1 = sin_N1_C.real();
662 T cos_N1 = -(N11+N22)/(2*R_N);
663 std::complex<T> sin_N2_C = -math::sqrt(sin_Nsqrt)/R_N;
664 T sin_N2 = sin_N2_C.real();
665 T cos_N2 = -(N11+N22)/(2*R_N);
666
667 T theta_0N = math::atan2(sin_N0,cos_N0);
668 T theta_1N = math::atan2(sin_N1,cos_N1);
669 T theta_2N = math::atan2(sin_N2,cos_N2);
670
671 T theta_1 = mod((theta_1N - theta_0N + theta_0E + pi),(2*pi));
672 T theta_2 = mod((theta_2N - theta_0N + theta_0E + pi),(2*pi));
673 if (math::isnan(theta_1) || math::isnan(theta_2) || R_N==1 || R_E == 1)
674 {
675 gsDebugVar(E);
676 gsDebugVar(C);
677 gsDebugVar(N);
678
679
680 gsDebugVar(R_E);
681
682 gsDebugVar(sin_Esqrt);
683
684 gsDebugVar(sin_E0);
685 gsDebugVar(sin_E1);
686 gsDebugVar(sin_E2);
687
688 gsDebugVar(cos_E0);
689 gsDebugVar(cos_E1);
690 gsDebugVar(cos_E2);
691
692 gsDebugVar(R_N);
693
694 gsDebugVar(sin_Nsqrt);
695
696 gsDebugVar(sin_N0);
697 gsDebugVar(sin_N1);
698 gsDebugVar(sin_N2);
699
700 gsDebugVar(cos_N0);
701 gsDebugVar(cos_N1);
702 gsDebugVar(cos_N2);
703
704
705 gsDebugVar(theta_1N);
706 gsDebugVar(theta_0N);
707 gsDebugVar(theta_0E);
708 }
709
710 T Q_E_min = theta_1E - theta_0E;
711 T Q_E_max = theta_2E - theta_0E;
712 T Q_S_min, Q_S_max;
713
714 T min, max;
715 if (theta_1 < theta_2)
716 {
717 Q_S_min = theta_1 - pi - theta_0E;
718 Q_S_max = theta_2 - pi - theta_0E;
719
720 // set difference:
721 min = math::max(Q_E_min,Q_S_min);
722 max = math::min(Q_E_max,Q_S_max);
723 min /=2;
724 max /=2;
725 // gsDebugVar("One interval");
726 }
727 else if (theta_1 > theta_2)
728 {
729 // try first interval
730 Q_S_min = theta_1 - pi - theta_0E;
731 Q_S_max = pi - theta_0E;
732
733 // set difference:
734 min = math::max(Q_E_min,Q_S_min);
735 max = math::min(Q_E_max,Q_S_max);
736 min /=2;
737 max /=2;
738 if (obj.eval(min)*obj.eval(max)>0)
739 {
740 // second interval
741 Q_S_min = - pi - theta_0E;
742 Q_S_max = theta_2 - pi - theta_0E;
743 // set difference:
744 min = math::max(Q_E_min,Q_S_min);
745 max = math::min(Q_E_max,Q_S_max);
746 min /=2;
747 max /=2;
748 // // if (obj.eval(min)*obj.eval(max)>0)
749 // // {
750 // // Q_S_min = theta_1 - pi - theta_0E;
751 // // Q_S_max = pi - theta_0E;
752 // // gsDebugVar(obj.eval(math::max(Q_E_min,Q_S_min)/2));
753 // // gsDebugVar(obj.eval(math::min(Q_E_max,Q_S_max)/2));
754 // // gsDebugVar(obj.eval(math::max(Q_E_min,Q_S_min)/2)*obj.eval(math::min(Q_E_max,Q_S_max)/2));
755 // // Q_S_min = - pi - theta_0E;
756 // // Q_S_max = theta_2 - pi - theta_0E;
757 // // gsDebugVar(obj.eval(math::max(Q_E_min,Q_S_min)/2));
758 // // gsDebugVar(obj.eval(math::min(Q_E_max,Q_S_max)/2));
759 // // gsDebugVar(obj.eval(math::max(Q_E_min,Q_S_min)/2)*obj.eval(math::min(Q_E_max,Q_S_max)/2));
760
761 // // GISMO_ERROR("Root is not found??");
762 // // }
763 }
764 // gsDebugVar("Two intervals");
765 }
766 else
767 {
768 min = 0;
769 max = 0;
770 // gsWarn<<"Interval not specified\n";
771 // GISMO_ERROR("Don't know what to do!");
772 }
773 result(0) = min;
774 result(1) = max;
775 return result;
776}
777
778template <short_t dim, class T, bool linear >
779gsMatrix<T> gsMaterialMatrixTFT<dim,T,linear>::_compute_TF(const index_t patch, const gsMatrix<T> & u, const gsMatrix<T> & z) const
780{
781 if (m_options.getSwitch("Explicit"))
782 {
783 function_ptr tmp = m_materialMat->getDeformed();
784 m_materialMat->setDeformed(m_defpatches0);
785 gsMatrix<T> TF = m_materialMat->eval3D_tensionfield(patch,u,z,MaterialOutput::TensionField);
786 m_materialMat->setDeformed(tmp);
787 return TF;
788 }
789 else
790 return m_materialMat->eval3D_tensionfield(patch,u,z,MaterialOutput::TensionField);
791}
792
793template <short_t dim, class T, bool linear >
794gsMatrix<T> gsMaterialMatrixTFT<dim,T,linear>::_compute_TF(const index_t patch, const gsVector<T>& u, const T & z) const
795{
796 gsMatrix<T> zmat(1,1);
797 zmat<<z;
798 return _compute_TF(patch,u,zmat);
799}
800
801
802template <short_t dim, class T, bool linear >
803gsMatrix<T> gsMaterialMatrixTFT<dim,T,linear>::_compute_E(const T theta, const gsMatrix<T> & C, const gsMatrix<T> & S, const gsMatrix<T> & E) const
804{
805 gsMatrix<T> result = E;
806
807 T n1 = math::cos(theta);
808 T n2 = math::sin(theta);
809 gsVector<T,3> n1_vec; n1_vec<<n1*n1, n2*n2, 2*n1*n2;
810
811 // Gamma
812 T gamma = - ( n1_vec.transpose() * S ).value() / ( n1_vec.transpose() * C * n1_vec ).value();
813
814 gsMatrix<T> Ew = gamma * n1_vec;
815 result += Ew;
816 return result;
817}
818
819template <short_t dim, class T, bool linear >
820gsMatrix<T> gsMaterialMatrixTFT<dim,T,linear>::_compute_S(const T theta, const gsMatrix<T> & C, const gsMatrix<T> & S) const
821{
822 gsMatrix<T> result = S;
823
824 T n1 = math::cos(theta);
825 T n2 = math::sin(theta);
826 gsVector<T,3> n1_vec; n1_vec<<n1*n1, n2*n2, 2*n1*n2;
827
828 // Gamma
829 T gamma = - ( n1_vec.transpose() * S ).value() / ( n1_vec.transpose() * C * n1_vec ).value();
830
831 gsMatrix<T> Ew = gamma * n1_vec;
832 result += C * Ew;
833 return result;
834}
835
836template <short_t dim, class T, bool linear >
837gsMatrix<T> gsMaterialMatrixTFT<dim,T,linear>::_compute_C(const T theta, const gsMatrix<T> & C, const gsMatrix<T> & S) const
838{
839 GISMO_ASSERT(C .rows()==3,"Number of rows of C must be 3, but is "<<C.rows());
840 GISMO_ASSERT(C .cols()==3,"Number of cols of C must be 3, but is "<<C.cols());
841 GISMO_ASSERT(S .rows()==3,"Number of rows of S must be 3, but is "<<S.rows());
842 GISMO_ASSERT(S .cols()==1,"Number of cols of S must be 1, but is "<<S.cols());
843
844 gsMatrix<T> result = C;
845
846 T n1 = math::cos(theta);
847 T n2 = math::sin(theta);
848 T m1 = -math::sin(theta);
849 T m2 = math::cos(theta);
850 gsVector<T,3> n1_vec; n1_vec<<n1*n1, n2*n2, 2*n1*n2;
851 gsVector<T,3> n2_vec; n2_vec<<m1*n1, m2*n2, m1*n2+m2*n1;
852 gsVector<T,3> n4_vec; n4_vec<<m1*m1, m2*m2, 2*m1*m2;
853
854 // Gamma
855 T gamma = - ( n1_vec.transpose() * S ).value() / ( n1_vec.transpose() * C * n1_vec ).value();
856
857 // dGamma/dE
858 gsMatrix<T> I(3,3); I.setIdentity();
859 gsMatrix<T> dgammadE = - ( n1_vec.transpose() * C) / ( n1_vec.transpose() * C * n1_vec ).value();
860 dgammadE.transposeInPlace();
861
862 // dGamma/dT
863 T tmp = ( n2_vec.transpose() * C * n1_vec ).value() / ( n1_vec.transpose() * C * n1_vec ).value();
864 T dgammadT = - 2 * gamma * tmp;
865
866 // dF/dE
867 // gsMatrix<T> dfdE = ( n2_vec.transpose() - tmp * n1_vec.transpose() ) * C;
868 gsMatrix<T> dfdE = n2_vec.transpose() * C - dgammadE.transpose() * (n2_vec.transpose() * C * n1_vec).value();
869 dfdE.transposeInPlace();
870
871 // dF/dT
872 // NAKASHINO
873 // T dfdT = (n4_vec.transpose() * ( S + C * gamma * n1_vec )).value() +
874 // 2 * gamma * ( n2_vec.transpose() * C * n2_vec - math::pow(( n1_vec.transpose() * C * n2_vec ).value(),2) / ( n1_vec.transpose() * C * n1_vec ).value() );
875 T dfdT = (n4_vec.transpose() * S).value() + dgammadT * (n2_vec.transpose() * C * n1_vec).value() +
876 gamma * ( (n4_vec.transpose() * C * n1_vec).value() + 2*(n2_vec.transpose() * C * n2_vec).value());
877
878
879 // dT/dE
880 gsMatrix<T> dTdE = - dfdE / dfdT;
881
882 result += C * ( n1_vec * dgammadE.transpose() + dgammadT * n1_vec * dTdE.transpose() + 2*gamma * n2_vec * dTdE.transpose() );
883 return result;
884}
885
886template <short_t dim, class T, bool linear >
887gsMatrix<T> gsMaterialMatrixTFT<dim,T,linear>::_compute_C(const T theta, const gsMatrix<T> & C, const gsMatrix<T> & S, const gsMatrix<T> & dC) const
888{
889 GISMO_ASSERT(C .rows()==3,"Number of rows of C must be 3, but is "<<C.rows());
890 GISMO_ASSERT(C .cols()==3,"Number of cols of C must be 3, but is "<<C.cols());
891 GISMO_ASSERT(S .rows()==3,"Number of rows of S must be 3, but is "<<S.rows());
892 GISMO_ASSERT(S .cols()==1,"Number of cols of S must be 1, but is "<<S.cols());
893 GISMO_ASSERT(dC.cols()==3,"Number of cols of dC must be 3, but is "<<dC.cols());
894 GISMO_ASSERT(dC.rows()==9,"Number of rows of dC must be 1, but is "<<dC.rows());
895
896 gsMatrix<T> result = C;
897
898 T n1 = math::cos(theta);
899 T n2 = math::sin(theta);
900 T m1 = -math::sin(theta);
901 T m2 = math::cos(theta);
902 gsVector<T,3> n1_vec; n1_vec<<n1*n1, n2*n2, 2*n1*n2;
903 gsVector<T,3> n2_vec; n2_vec<<m1*n1, m2*n2, m1*n2+m2*n1;
904 gsVector<T,3> n4_vec; n4_vec<<m1*m1, m2*m2, 2*m1*m2;
905
906 // Derivative of C to C, dC/dC, times N1
907 gsMatrix<T> dCN1(3,3);
908 for (index_t d = 0; d!=dC.cols(); d++)
909 {
910 dCN1.col(d) = dC.reshapeCol(d,3,3) * n1_vec;
911 }
912
913 // Gamma
914 T gamma = - ( n1_vec.transpose() * S ).value() / ( n1_vec.transpose() * C * n1_vec ).value();
915
916 // dGamma/dE
917 gsMatrix<T> I(3,3); I.setIdentity();
918 gsMatrix<T> dgammadE = - ( n1_vec.transpose() * C + gamma * ( n1_vec.transpose() * dCN1) ) / ( n1_vec.transpose() * C * n1_vec ).value();
919 // dgammadE.transposeInPlace();
920
921 // dGamma/dT
922 T tmp = ( n2_vec.transpose() * C * n1_vec ).value() / ( n1_vec.transpose() * C * n1_vec ).value();
923 T dgammadT = - 2 * gamma * tmp;
924
925 // dF/dE
926 // gsMatrix<T> dfdE = ( n2_vec.transpose() - tmp * n1_vec.transpose() ) * C;
927 gsMatrix<T> dfdE = n2_vec.transpose() * C + dgammadE * (n2_vec.transpose() * C * n1_vec).value() + gamma * ( n2_vec.transpose() * dCN1);
928 // dfdE.transposeInPlace();
929
930 // dF/dT
931 // NAKASHINO
932 // T dfdT = (n4_vec.transpose() * ( S + C * gamma * n1_vec )).value() +
933 // 2 * gamma * ( n2_vec.transpose() * C * n2_vec - math::pow(( n1_vec.transpose() * C * n2_vec ).value(),2) / ( n1_vec.transpose() * C * n1_vec ).value() );
934 T dfdT = (n4_vec.transpose() * S).value() + dgammadT * (n2_vec.transpose() * C * n1_vec).value() +
935 gamma * ( (n4_vec.transpose() * C * n1_vec).value() + 2*(n2_vec.transpose() * C * n2_vec).value());
936
937 // dT/dE
938 gsMatrix<T> dTdE = - dfdE / dfdT;
939
940 // d(gamma*n1)dE
941 gsMatrix<T> dgamman1dE = ( n1_vec * dgammadE + dgammadT * n1_vec * dTdE + 2*gamma * n2_vec * dTdE );
942 result += dgamman1dE.transpose() * C;
943 result += gamma * dCN1;
944 return result;
945}
946} // end namespace
947
948using namespace gismo;
949
950template <typename T>
951T mod(const T x, const T y)
952{
953 return x - math::floor(x/y)*y;
954}
955
956template <typename T>
957class objective : public gsFunction<T>
958{
959public:
960 objective(const gsMatrix<T> & C, const gsMatrix<T> & S)
961 :
962 m_C(C),
963 m_S(S)
964 {
965 m_support.resize(1,2);
966 m_support<<0,2*3.1415926535;
967 }
968
969 void eval_into(const gsMatrix<T>& u, gsMatrix<T>& result) const
970 {
971 result.resize(1,u.cols());
972 gsVector<T,3> n1_vec, n2_vec;
973 T theta, gamma;
974 T n1, n2, m1, m2;
975
976 for (index_t k = 0; k!=u.cols(); k++)
977 {
978 theta = u(0,k);
979 n1 = math::cos(theta);
980 n2 = math::sin(theta);
981 m1 = -math::sin(theta);
982 m2 = math::cos(theta);
983
984 n1_vec<<n1*n1, n2*n2, 2*n1*n2;
985 n2_vec<<m1*n1, m2*n2, m1*n2+m2*n1;
986
987 gamma = - ( n1_vec.transpose() * m_S ).value() / ( n1_vec.transpose() * m_C * n1_vec ).value();
988
989 T res = (n2_vec.transpose() * m_S).value() + gamma * (n2_vec.transpose() * m_C * n1_vec).value();
990
991 result(0,k) = res;
992 }
993 }
994
995 T eval(const T & theta) const
996 {
997 gsMatrix<T> t(1,1);
998 t<<theta;
999 gsMatrix<T> res;
1000 this->eval_into(t,res);
1001 return res(0,0);
1002 }
1003
1004 short_t domainDim() const
1005 {
1006 return 1;
1007 }
1008
1009 short_t targetDim() const
1010 {
1011 return 1;
1012 }
1013
1014 gsMatrix<T> support() const
1015 {
1016 return m_support;
1017 }
1018
1019 void deriv_into2(const gsMatrix<T>& u, gsMatrix<T>& result) const
1020 {
1021 result.resize(1,u.cols());
1022 gsVector<T,3> n1_vec, n2_vec, n3_vec, n4_vec;
1023 T theta, gamma, dgammadT;
1024 T n1, n2, m1, m2;
1025
1026 for (index_t k = 0; k!=u.cols(); k++)
1027 {
1028 theta = u(0,k);
1029 n1 = math::cos(theta);
1030 n2 = math::sin(theta);
1031 m1 = -math::sin(theta);
1032 m2 = math::cos(theta);
1033
1034 n1_vec<<n1*n1, n2*n2, 2*n1*n2;
1035 n2_vec<<m1*n1, m2*n2, m1*n2+m2*n1;
1036 n3_vec<<m1*m1-n1*n1, m2*m2-n2*n2, 2*(m1*m2-n1*n2);
1037 n4_vec<<m1*m1, m2*m2, 2*m1*m2;
1038
1039 gamma = - ( n1_vec.transpose() * m_S ).value() / ( n1_vec.transpose() * m_C * n1_vec ).value();
1040 T tmp = ( n2_vec.transpose() * m_C * n1_vec ).value() / ( n1_vec.transpose() * m_C * n1_vec ).value();
1041 dgammadT = - 2 * gamma * tmp;
1042 // gsDebugVar(dgammadT);
1043
1044 T res = (n3_vec.transpose() * m_S).value() + dgammadT * (n2_vec.transpose() * m_C * n1_vec).value()
1045 + gamma * (n3_vec.transpose() * m_C * n1_vec).value() + 2*gamma * (n2_vec.transpose() * m_C * n2_vec).value();
1046 // T res = (n2_vec.transpose() * m_S).value() + gamma * (n2_vec.transpose() * m_C * n1_vec).value();
1047
1048 result(0,k) = res;
1049 }
1050 }
1051
1052 void setSupport(const gsMatrix<T> supp)
1053 {
1054 m_support = supp;
1055 }
1056
1069 bool findRootBisection(T & /*f*/, const T & A, const T & B, T & x, const T & t = 1e-12, const index_t & itmax = 1000)
1070 {
1071 T a = A;
1072 T b = B;
1073 T c = (b-a)/2;
1074 index_t it = 0;
1075 while (math::abs(b-a) > t && it++ < itmax)
1076 {
1077 if (this->eval(c)==0.0)
1078 break;
1079 else if (this->eval(c)*this->eval(a) < 0)
1080 b = c;
1081 else
1082 a = c;
1083 }
1084 x = c;
1085 if (math::abs(b-a) < t)
1086 return true;
1087 else
1088 return false;
1089 }
1090
1105 bool findRootBrent(T & f, const T & a, const T & b, T & x, const T & t = 1e-12, const index_t & itmax = 1000)
1106 {
1107 /*
1108 * Implementation of Brent's method for finding the root \a c of a scalar equation \a *this within the interval [a,b]
1109 *
1110 * Adopted from: https://people.math.sc.edu/Burkardt/cpp_src/brent/brent.html
1111 * The brent.cpp file, function zero( )
1112 * https://people.math.sc.edu/Burkardt/cpp_src/brent/brent.cpp
1113 */
1114 T c;
1115 T d;
1116 T e;
1117 T fa;
1118 T fb;
1119 T fc;
1120 T m;
1121 T macheps;
1122 T p;
1123 T q;
1124 T r;
1125 T s;
1126 T sa;
1127 T sb;
1128 T tol;
1129//
1130// Make local copies of A and B.
1131//
1132 sa = a;
1133 sb = b;
1134 fa = this->eval( sa );
1135 fb = this->eval( sb );
1136
1137 if (fa*fb >=0)
1138 return false;
1139 else if (fa > fb)
1140 std::swap(fa,fb);
1141
1142 c = sa;
1143 fc = fa;
1144 e = sb - sa;
1145 d = e;
1146
1147 macheps = std::numeric_limits<T>::epsilon();
1148
1149 for ( index_t it = 0; it!=itmax; it++ )
1150 {
1151 if ( std::fabs ( fc ) < std::fabs ( fb ) )
1152 {
1153 sa = sb;
1154 sb = c;
1155 c = sa;
1156 fa = fb;
1157 fb = fc;
1158 fc = fa;
1159 }
1160
1161 tol = 2.0 * macheps * std::fabs ( sb ) + t;
1162 m = 0.5 * ( c - sb );
1163
1164 if ( std::fabs ( m ) <= tol || fb == 0.0 )
1165 {
1166 x = sb;
1167 f = std::fabs ( m );
1168 return true;
1169 // break;
1170 }
1171
1172 if ( std::fabs ( e ) < tol || std::fabs ( fa ) <= std::fabs ( fb ) )
1173 {
1174 e = m;
1175 d = e;
1176 }
1177 else
1178 {
1179 s = fb / fa;
1180
1181 if ( sa == c )
1182 {
1183 p = 2.0 * m * s;
1184 q = 1.0 - s;
1185 }
1186 else
1187 {
1188 q = fa / fc;
1189 r = fb / fc;
1190 p = s * ( 2.0 * m * q * ( q - r ) - ( sb - sa ) * ( r - 1.0 ) );
1191 q = ( q - 1.0 ) * ( r - 1.0 ) * ( s - 1.0 );
1192 }
1193
1194 if ( 0.0 < p )
1195 {
1196 q = - q;
1197 }
1198 else
1199 {
1200 p = - p;
1201 }
1202
1203 s = e;
1204 e = d;
1205
1206 if ( 2.0 * p < 3.0 * m * q - std::fabs ( tol * q ) &&
1207 p < std::fabs ( 0.5 * s * q ) )
1208 {
1209 d = p / q;
1210 }
1211 else
1212 {
1213 e = m;
1214 d = e;
1215 }
1216 }
1217 sa = sb;
1218 fa = fb;
1219
1220 if ( tol < std::fabs ( d ) )
1221 {
1222 sb = sb + d;
1223 }
1224 else if ( 0.0 < m )
1225 {
1226 sb = sb + tol;
1227 }
1228 else
1229 {
1230 sb = sb - tol;
1231 }
1232
1233 fb = this->eval( sb );
1234
1235 if ( ( 0.0 < fb && 0.0 < fc ) || ( fb <= 0.0 && fc <= 0.0 ) )
1236 {
1237 c = sa;
1238 fc = fa;
1239 e = sb - sa;
1240 d = e;
1241 }
1242 }
1243
1244
1245 GISMO_ERROR("Brent's method did not find a root");
1246 return false;
1247 }
1248
1249private:
1250 const gsMatrix<T> m_C;
1251 const gsMatrix<T> m_S;
1252 mutable gsMatrix<T> m_support;
1253};
gismo::gsAsMatrix
Creates a mapped object or data pointer to a matrix without copying data.
Definition gsAsMatrix.h:32
gismo::gsFunctionSet::eval
gsMatrix< T > eval(const gsMatrix< T > &u) const
Evaluate the function,.
Definition gsFunctionSet.hpp:120
gismo::gsFunction
A function from a n-dimensional domain to an m-dimensional image.
Definition gsFunction.h:60
gismo::gsFunction::domainDim
virtual short_t domainDim() const=0
Dimension of the (source) domain.
gismo::gsFunction::targetDim
virtual short_t targetDim() const
Dimension of the target space.
Definition gsFunctionSet.h:595
gismo::gsFunction::eval_into
virtual void eval_into(const gsMatrix< T > &u, gsMatrix< T > &result) const =0
Evaluate the function at points u into result.
gismo::gsMaterialMatrixTFT
This class defines a linear material.
Definition gsMaterialMatrixTFT.h:39
gismo::gsMaterialMatrixTFT::eval3D_vector
gsMatrix< T > eval3D_vector(const index_t patch, const gsMatrix< T > &u, const gsMatrix< T > &z, enum MaterialOutput out=MaterialOutput::Generic) const override
See gsMaterialMatrixBase for details.
Definition gsMaterialMatrixTFT.hpp:55
gismo::gsMaterialMatrixTFT::eval3D_strain
gsMatrix< T > eval3D_strain(const index_t patch, const gsMatrix< T > &u, const gsMatrix< T > &z) const override
See gsMaterialMatrixBase for details.
Definition gsMaterialMatrixTFT.hpp:117
gismo::gsMaterialMatrixTFT::eval3D_CauchyStress
gsMatrix< T > eval3D_CauchyStress(const index_t patch, const gsMatrix< T > &u, const gsMatrix< T > &z, enum MaterialOutput out) const override
See gsMaterialMatrixBase for details.
Definition gsMaterialMatrixTFT.hpp:199
gismo::gsMaterialMatrixTFT::eval3D_gamma
gsMatrix< T > eval3D_gamma(const index_t patch, const gsMatrix< T > &u, const gsMatrix< T > &z, enum MaterialOutput out=MaterialOutput::Generic) const override
See gsMaterialMatrixBase for details.
Definition gsMaterialMatrixTFT.hpp:336
gismo::gsMaterialMatrixTFT::eval3D_pstrain
gsMatrix< T > eval3D_pstrain(const index_t patch, const gsMatrix< T > &u, const gsMatrix< T > &z) const override
See gsMaterialMatrixBase for details.
Definition gsMaterialMatrixTFT.hpp:89
gismo::gsMaterialMatrixTFT::print
std::ostream & print(std::ostream &os) const override
See gsMaterialMatrixBase for details.
Definition gsMaterialMatrixTFT.hpp:29
gismo::gsMaterialMatrixTFT::eval3D_theta
gsMatrix< T > eval3D_theta(const index_t patch, const gsMatrix< T > &u, const gsMatrix< T > &z, enum MaterialOutput out=MaterialOutput::Generic) const override
See gsMaterialMatrixBase for details.
Definition gsMaterialMatrixTFT.hpp:304
gismo::gsMaterialMatrixTFT::eval3D_tensionfield
gsMatrix< T > eval3D_tensionfield(const index_t patch, const gsMatrix< T > &u, const gsMatrix< T > &z, enum MaterialOutput out=MaterialOutput::Generic) const override
See gsMaterialMatrixBase for details.
Definition gsMaterialMatrixTFT.hpp:298
gismo::gsMaterialMatrixTFT::eval_theta
gsMatrix< T > eval_theta(const gsMatrix< T > &Cs, const gsMatrix< T > &Ns, const gsMatrix< T > &Es) const
Computes theta.
Definition gsMaterialMatrixTFT.hpp:455
gismo::gsMaterialMatrixTFT::_eval3D_matrix_impl
std::enable_if< _linear, gsMatrix< T > >::type _eval3D_matrix_impl(const index_t patch, const gsMatrix< T > &u, const gsMatrix< T > &z, enum MaterialOutput out=MaterialOutput::Generic) const
Provides an implementation of eval3D_matrix for gsMaterialMatrixLinear.
Definition gsMaterialMatrixTFT.hpp:371
gismo::gsMaterialMatrixTFT::eval3D_stress
gsMatrix< T > eval3D_stress(const index_t patch, const gsMatrix< T > &u, const gsMatrix< T > &z, enum MaterialOutput out) const override
See gsMaterialMatrixBase for details.
Definition gsMaterialMatrixTFT.hpp:155
gismo::gsMaterialMatrixTFT::eval3D_pstress
gsMatrix< T > eval3D_pstress(const index_t patch, const gsMatrix< T > &u, const gsMatrix< T > &z, enum MaterialOutput out=MaterialOutput::Generic) const override
Definition gsMaterialMatrixTFT.hpp:61
gismo::gsMaterialMatrixTFT::eval3D_matrix
gsMatrix< T > eval3D_matrix(const index_t patch, const gsMatrix< T > &u, const gsMatrix< T > &z, enum MaterialOutput out=MaterialOutput::Generic) const override
See gsMaterialMatrixBase for details.
Definition gsMaterialMatrixTFT.hpp:42
gismo::gsMatrix
A matrix with arbitrary coefficient type and fixed or dynamic size.
Definition gsMatrix.h:41
gismo::gsMatrix::reshapeCol
gsAsMatrix< T, Dynamic, Dynamic > reshapeCol(index_t c, index_t n, index_t m)
Returns column c of the matrix resized to n x m matrix This function assumes that the matrix is size ...
Definition gsMatrix.h:231
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::gsVector
A vector with arbitrary coefficient type and fixed or dynamic size.
Definition gsVector.h:37
gismo::MaterialOutput
MaterialOutput
This class describes the output type.
Definition gsMaterialMatrixUtils.h:99
short_t
#define short_t
Definition gsConfig.h:35
index_t
#define index_t
Definition gsConfig.h:32
GISMO_ERROR
#define GISMO_ERROR(message)
Definition gsDebug.h:118
gsWarn
#define gsWarn
Definition gsDebug.h:50
GISMO_ASSERT
#define GISMO_ASSERT(cond, message)
Definition gsDebug.h:89
gsFunction.h
Provides declaration of Function abstract interface.
gsMaterialMatrixUtils.h
Provides material matrix utilities.
gismo
The G+Smo namespace, containing all definitions for the library.