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