G+Smo  25.01.0
Geometry + Simulation Modules
 
Loading...
Searching...
No Matches
gsALMCrisfield.hpp
Go to the documentation of this file.
1
14#pragma once
15
16#include <typeinfo>
17#include <gsStructuralAnalysis/src/gsALMSolvers/gsALMHelper.h>
18
19namespace gismo
20{
21
22template <class T>
24{
25 Base::defaultOptions();
26 m_options.addReal("Scaling","Set Scaling factor Phi",-1);
27 m_options.addInt ("AngleMethod","Angle determination method: 0 = Previous step; 1 = Previous iteration",angmethod::Step);
28}
29
30template <class T>
32{
33 Base::getOptions();
34 m_phi = m_options.getReal("Scaling");
35 m_phi_user = m_phi == -1 ? false : true;
36
37 m_angleDetermine = m_options.getInt ("AngleMethod");
38}
39
40template <class T>
42{
43 m_numDof = m_forcing.size();
44 m_DeltaU = m_U = gsVector<T>::Zero(m_numDof);
45 m_DeltaL = m_L = 0.0;
46
47 m_DeltaUold = gsVector<T>::Zero(m_numDof);
48 m_DeltaLold = 0.0;
49}
50
51// ------------------------------------------------------------------------------------------------------------
52// ---------------------------------------Crisfield's method---------------------------------------------------
53// ------------------------------------------------------------------------------------------------------------
54
55template <class T>
57{
58 m_jacMat = computeJacobian();
59 this->factorizeMatrix(m_jacMat);
60 computeUt(); // rhs does not depend on solution
61 computeUbar(); // rhs contains residual and should be computed every time
62
63}
64
65template <class T>
67{
68 m_jacMat = computeJacobian();
69 this->factorizeMatrix(m_jacMat);
70 computeUt(); // rhs does not depend on solution
71}
72
73template <class T>
75{
76 computeUbar(); // rhs contains residual and should be computed every time
77
78 // Compute next solution
79 m_eta = 1.0;
80
81 T lamold = m_deltaL;
82 computeLambdas();
83
84 // Relaxation against oscillating load factor
85 if (( (lamold*m_deltaL < 0) && (abs(m_deltaL) <= abs(lamold) ) ) && m_relax != 1.0 )
86 {
87 m_note += "\t relaxated solution!";
88 m_deltaU = m_relax * (m_deltaL*m_deltaUt + m_eta*m_deltaUbar);
89 m_deltaL = m_relax * m_deltaL;
90 }
91
92 m_DeltaU += m_deltaU;
93 m_DeltaL += m_deltaL;
94
95 if (m_angleDetermine == angmethod::Iteration)
96 {
97 m_DeltaUold = m_DeltaU;
98 m_DeltaLold = m_DeltaL;
99 }
100}
101
102template <class T>
104{
105 m_DeltaU = m_deltaUbar = m_deltaUt = gsVector<T>::Zero(m_numDof);
106 m_DeltaL = m_deltaL = 0.0;
107 m_eta = 1.0;
108}
109
110template <class T>
112{
113 m_jacMat = computeJacobian();
114 this->factorizeMatrix(m_jacMat);
115 m_deltaUt = this->solveSystem(m_forcing);
116
117 // Choose Solution
118 if (m_DeltaUold.dot(m_DeltaUold) == 0 && m_DeltaLold*m_DeltaLold == 0) // no information about previous step.
119 {
120 // gsWarn<<"Different predictor!!!!!\n";
121 m_note+= "predictor\t";
122 m_deltaL = m_arcLength / math::pow(2*( m_deltaUt.dot(m_deltaUt) ) , 0.5);
123 m_deltaU = m_deltaUbar + m_deltaL*m_deltaUt;
124 if (!m_phi_user)
125 m_phi = math::pow( m_deltaUt.dot(m_deltaUt) / m_forcing.dot(m_forcing),0.5);
126 m_note += " phi=" + std::to_string(m_phi);
127
128 }
129 // {
130 // m_note+= "predictor\t";
131
132 // if (!m_phi_user)
133 // m_phi = math::pow( m_deltaUt.dot(m_deltaUt) / m_forcing.dot(m_forcing),0.5);
134
135 // // m_deltaL = m_arcLength * DL / math::sqrt( 2*( m_deltaUt.dot(m_deltaUt) + m_DeltaL*DL ) );
136 // m_deltaL = m_arcLength / math::sqrt( ( m_deltaUt.dot(m_deltaUt) + m_phi*m_phi ) );
137
138 // // m_deltaU = m_arcLength * m_deltaUt / math::sqrt( m_deltaUt.dot(m_deltaUt) + m_DeltaL*DL );
139 // m_deltaU = m_deltaL*m_deltaUt;
140
141
142 // m_note += " phi=" + std::to_string(m_phi);
143
144 // // m_DeltaUold = m_deltaU;
145 // // m_DeltaLold = m_deltaL;
146 // }
147 else // previous point is not in the origin
148 {
149 if (!m_phi_user)
150 m_phi = math::pow(m_U.dot(m_U)/( math::pow(m_L,2) * m_forcing.dot(m_forcing) ),0.5);
151 m_note += " phi=" + std::to_string(m_phi);
152 computeLambdaMU();
153 }
154
155 // Compute Temporary updates of DeltaL and DeltaU
156 m_DeltaU += m_deltaU;
157 m_DeltaL += m_deltaL;
158
159 if (m_angleDetermine == angmethod::Iteration || m_angleDetermine == angmethod::Predictor)
160 {
161 m_DeltaUold = m_DeltaU;
162 m_DeltaLold = m_DeltaL;
163 }
164}
165
166template <class T>
168{
169 GISMO_ASSERT(m_Uguess.rows()!=0 && m_Uguess.cols()!=0,"Guess is empty");
170
171 m_jacMat = computeJacobian();
172 this->factorizeMatrix(m_jacMat);
173 m_deltaUt = this->solveSystem(m_forcing);
174 if (!m_phi_user)
175 m_phi = math::pow( m_deltaUt.dot(m_deltaUt) / m_forcing.dot(m_forcing),0.5);
176 m_note += " phi=" + std::to_string(m_phi);
177
178 //
179 m_DeltaUold = -(m_Uguess - m_U);
180 m_DeltaLold = -(m_Lguess - m_L);
181
182 // m_DeltaUold *= m_arcLength / math::sqrt( m_deltaU.dot(m_deltaU));
183 // m_DeltaLold *= m_arcLength / math::sqrt( m_deltaU.dot(m_deltaU));
184
185 computeLambdaMU();
186
187 m_DeltaU = m_deltaU;
188 m_DeltaL = m_deltaL;
189
190 m_Uguess.resize(0);
191}
192
193template <class T>
195{
196 m_converged = true;
197 m_Uprev = m_U;
198 m_Lprev = m_L;
199 m_U += m_DeltaU;
200 m_L += m_DeltaL;
201 if (m_angleDetermine == angmethod::Step)
202 {
203 m_DeltaUold = m_DeltaU;
204 m_DeltaLold = m_DeltaL;
205 }
206}
207
208// ------------------------------------------------------------------------------------------------------------
209// ---------------------------------------Lambda computations--------------------------------------------------
210// ------------------------------------------------------------------------------------------------------------
211
212template <class T>
213void gsALMCrisfield<T>::computeLambdasSimple() //Ritto-Corrêa et al. 2008
214{
215 T A0 = math::pow(m_phi,2)* m_forcing.dot(m_forcing); // see Lam et al. 1991,
216
217 m_a0 = m_deltaUt.dot(m_deltaUt) + A0;
218 m_b0 = 2*( m_deltaUt.dot(m_DeltaU) + m_DeltaL * A0 );
219 m_b1 = 2*( m_deltaUbar.dot(m_deltaUt) );
220 m_c0 = m_DeltaU.dot(m_DeltaU) + m_DeltaL*m_DeltaL * A0 - math::pow(m_arcLength,2);
221 m_c1 = 2*( m_DeltaU.dot(m_deltaUbar) );
222 m_c2 = m_deltaUbar.dot(m_deltaUbar);
223
225 m_alpha1 = m_a0;
226 m_alpha2 = m_b0 + m_eta*m_b1;
227 m_alpha3 = m_c0 + m_eta*m_c1 + m_eta*m_eta*m_c2;
228
229 m_discriminant = math::pow(m_alpha2 ,2) - 4 * m_alpha1 * m_alpha3;
230 // m_note += "\t D = " + std::to_string(m_discriminant);
231}
232
233template <class T>
235{
236 m_alpha1 = m_a0;
237 m_alpha2 = m_b0 + m_eta*m_b1;
238 m_alpha3 = m_c0 + m_eta*m_c1 + m_eta*m_eta*m_c2;
239
240 // Lam 1992
241 // m_discriminant = math::pow(m_alpha2 ,2) - 4 * m_alpha1 * m_alpha3;
242 // m_deltaLs[0] = (-m_alpha2 + math::sqrt(m_discriminant))/(2*m_alpha1);
243 // m_deltaLs[1] = (-m_alpha2 - math::sqrt(m_discriminant))/(2*m_alpha1);
244
245 // Zhou 1995
246 m_deltaLs[0] = (-m_alpha2 )/(2*m_alpha1);
247 m_deltaLs[1] = (-m_alpha2 )/(2*m_alpha1);
248}
249
250template <class T>
252{
253 m_alpha1 = m_b1*m_b1 - 4.0*m_a0*m_c2;
254 m_alpha2 = 2.0*m_b0*m_b1 - 4.0*m_a0*m_c1;
255 m_alpha3 = m_b0*m_b0 - 4.0*m_a0*m_c0;
256
257 m_discriminant = math::pow(m_alpha2 ,2.0) - 4.0 * m_alpha1 * m_alpha3;
258
259 gsVector<T> etas(2);
260 etas.setZero();
261 if (m_discriminant >= 0)
262 {
263 etas[0] = (-m_alpha2 + math::pow(m_discriminant,0.5))/(2.0*m_alpha1);
264 etas[1] = (-m_alpha2 - math::pow(m_discriminant,0.5))/(2.0*m_alpha1);
265
266 T eta1 = std::min(etas[0],etas[1]);
267 T eta2 = std::max(etas[0],etas[1]);
268 if (m_verbose) {gsInfo<<"eta 1 = "<<eta1<<"\t eta2 = "<<eta2<<"\n";}
269
270 // Approach of Zhou 1995
271 // m_eta = std::min(1.0,eta2);
272 // if (m_eta <= 0)
273 // gsInfo<<"Warning: both etas are non-positive!\n";
274 // if (m_eta <=0.5)
275 // {
276 // gsInfo<<"Warning: eta is small; bisecting step length\n";
277 // m_arcLength=m_arcLength/2.;
278 // }
279
280 // Approach of Lam 1992
281 T xi = 0.05*abs(eta2-eta1);
282 if (eta2<1.0)
283 m_eta = eta2-xi;
284 else if ( (eta2 > 1.0) && (-m_alpha2/m_alpha1 < 1.0) )
285 m_eta = eta2+xi;
286 else if ( (eta1 < 1.0) && (-m_alpha2/m_alpha1 > 1.0) )
287 m_eta = eta1-xi;
288 else if ( eta1 > 1.0 )
289 m_eta = eta1 + xi;
290
291 if (eta2<1.0)
292 m_note = m_note + " option 1";
293 else if ( (eta2 > 1.0) && (-m_alpha2/m_alpha1 < 1.0) )
294 m_note = m_note + " option 2";
295 else if ( (eta1 < 1.0) && (-m_alpha2/m_alpha1 > 1.0) )
296 m_note = m_note + " option 3";
297 else if ( eta1 > 1.0 )
298 m_note = m_note + " option 4";
299
300 // m_eta = eta2;
301 }
302 else
303 {
304 gsInfo<<"Discriminant was negative in modified method\n";
305 }
306}
307template <class T>
309{
310 T A0 = math::pow(m_phi,2)* m_forcing.dot(m_forcing); // see Lam et al. 1991
311 gsVector<T> DeltaUcr = m_DeltaU + m_deltaUbar;
312
313 // Compute internal loads from residual and
314 // gsVector<T> R = m_residualFun(m_U + DeltaUcr, m_L + m_DeltaL , m_forcing);
315 // gsVector<T> Fint = R + (m_L + m_DeltaL) * m_forcing;
316 gsVector<T> Fint = m_jacMat*(m_U+m_DeltaU);
317 T Lcr = Fint.dot(m_forcing)/m_forcing.dot(m_forcing);
318 T DeltaLcr = Lcr - m_L;
319
320 T arcLength_cr = math::pow( DeltaUcr.dot(DeltaUcr) + A0 * math::pow(DeltaLcr, 2.0) ,0.5);
321 T mu = m_arcLength/arcLength_cr;
322
323 m_deltaL = mu*DeltaLcr - m_DeltaL;
324 m_deltaU = mu*DeltaUcr - m_DeltaU;
325}
326
327template <class T>
329{
330 m_deltaLs.setZero(2);
331 computeLambdasSimple();
332 if (m_discriminant >= 0)
333 {
334 m_eta = 1.0;
335 m_deltaLs[0] = (-m_alpha2 + math::pow(m_discriminant,0.5))/(2*m_alpha1);
336 m_deltaLs[1] = (-m_alpha2 - math::pow(m_discriminant,0.5))/(2*m_alpha1);
337 computeLambdaDOT();
338 }
339 else
340 {
341 m_note += "\tC";
342 // Compute eta
343 computeLambdasModified();
344 gsDebugVar(m_discriminant);
345 if ((m_discriminant >= 0) && m_eta > 0.05)
346 {
347 // recompute lambdas with new eta
348 computeLambdasEta();
349 // gsInfo<<"2: dL1 = "<<m_deltaLs[0]<<"\tdL2 = "<<m_deltaLs[1]<<"\t eta = "<<m_eta<<"\n";
350 computeLambdaDOT();
351 // gsInfo<<"2: dL1 = "<<m_deltaL<<"\t m_deltaU.norm = "<<m_deltaU.norm()<<"\t eta = "<<m_eta<<"\n";
352 if (m_verbose) {gsInfo<<"Modified Complex Root Solve\n";}
353 }
354 else
355 {
356 // if the roots of the modified method are still complex, we use the following function (see Lam 1992, eq 13-17)
357 m_eta = 1.0;
358 computeLambdasComplex();
359 // gsInfo<<"3: dL1 = "<<m_deltaL<<"\t m_deltaU.norm = "<<m_deltaU.norm()<<"\t eta = "<<m_eta<<"\n";
360 if (m_verbose) {gsInfo<<"Simplified Complex Root Solve\n";}
361 // Note: no selection of roots is needed
362 }
363 }
364}
365
366template <class T>
368{
369 computeLambdas();
370
371 if (sign(m_DeltaL + m_deltaLs[0]) == sign(m_detKT))
372 m_deltaL = m_deltaLs[0];
373 else
374 m_deltaL = m_deltaLs[1];
375
376 // Compute update of U (NOTE: m_eta=1.0)
377 m_deltaU = m_deltaUbar + m_deltaL*m_deltaUt;
378
379 // gsInfo<<"\t\t Choice based on DETERMINANT. Options:\n";
380 // gsInfo<<"\t\t DeltaL = "<<m_DeltaL+m_deltaLs[0]<<" DeltaU.norm = "<<(m_DeltaU + m_deltaUbar + m_deltaLs[0]*m_deltaUt).norm()<<"\n";
381 // gsInfo<<"\t\t DeltaL = "<<m_DeltaL+m_deltaLs[1]<<" DeltaU.norm = "<<(m_DeltaU + m_deltaUbar + m_deltaLs[1]*m_deltaUt).norm()<<"\n";
382}
383
384template <class T>
386{
387 T A0 = math::pow(m_phi,2)* m_forcing.dot(m_forcing); // see Lam et al. 1991
388 index_t dir = sign(m_DeltaUold.dot(m_deltaUt) + A0*m_DeltaLold); // Feng et al. 1995 with H = \Psi^2
389 T denum = ( math::pow( m_deltaUt.dot(m_deltaUt) + A0 ,0.5) ); // Feng et al. 1995 with H = \Psi^2
390
391 T mu;
392 if (denum==0)
393 mu = m_arcLength;
394 else
395 mu = m_arcLength / denum;
396
397 m_deltaL = dir*mu;
398 m_deltaU = m_deltaL*m_deltaUt;
399}
400
401template <class T>
403{
404 gsVector<T> deltaU1, deltaU2;
405 deltaU1 = m_eta*m_deltaUbar + m_deltaUt*m_deltaLs[0];
406 deltaU2 = m_eta*m_deltaUbar + m_deltaUt*m_deltaLs[1];
407
408 // ---------------------------------------------------------------------------------
409 // Method by Ritto-Corea et al. 2008
410 T DOT1,DOT2;
411 DOT1 = m_deltaLs[0]*(m_DeltaUold.dot(m_deltaUt) + math::pow(m_phi,2)*m_DeltaLold);
412 DOT2 = m_deltaLs[1]*(m_DeltaUold.dot(m_deltaUt) + math::pow(m_phi,2)*m_DeltaLold);
413
414 if (DOT1 > DOT2)
415 {
416 m_deltaL = m_deltaLs[0];
417 m_deltaU = deltaU1;
418 }
419 else if (DOT1 < DOT2)
420 {
421 m_deltaL = m_deltaLs[1];
422 m_deltaU = deltaU2;
423 }
424 else
425 {
426 m_deltaL = m_deltaLs[0];
427 m_deltaU = deltaU1;
428 }
429
430 // ---------------------------------------------------------------------------------
431 // // Method by Crisfield 1981
432 // T DOT1,DOT2;
433 // DOT1 = (m_DeltaUold+deltaU1).dot(m_DeltaUold);
434 // DOT2 = (m_DeltaUold+deltaU2).dot(m_DeltaUold);
435
436 // DOT1 = (m_DeltaU+deltaU1).dot(m_DeltaUold);
437 // DOT2 = (m_DeltaU+deltaU2).dot(m_DeltaUold);
438
439 // if ((DOT1 > DOT2) && DOT2 <= 0)
440 // {
441 // m_deltaL = m_deltaLs[0];
442 // m_deltaU = deltaU1;
443 // }
444 // else if ((DOT1 < DOT2) && DOT1 <= 0)
445 // {
446 // m_deltaL = m_deltaLs[1];
447 // m_deltaU = deltaU2;
448 // }
449 // else if ((DOT1 >=0) && (DOT2 >=0))
450 // {
451 // T linsol = -m_alpha3/m_alpha2;
452 // T diff1 = abs(m_deltaLs[0]-linsol);
453 // T diff2 = abs(m_deltaLs[1]-linsol);
454 // m_note += "\t linear solution!";
455 // // m_note += "Linsol\t" + std::to_string(diff1) + "\t" + std::to_string(diff2) + "\t" + std::to_string(linsol) + "\t" + std::to_string(m_deltaLs[0]) + "\t" + std::to_string(m_deltaLs[1]) + "\n";
456 // if (diff1 > diff2)
457 // {
458 // m_deltaL = m_deltaLs[1];
459 // m_deltaU = deltaU2;
460 // }
461 // else
462 // {
463 // m_deltaL = m_deltaLs[0];
464 // m_deltaU = deltaU1;
465 // }
466 // }
467 // else
468 // {
469 // m_deltaL = m_deltaLs[0];
470 // m_deltaU = deltaU1;
471 // }
472
473}
474
475// ------------------------------------------------------------------------------------------------------------
476// ---------------------------------------Output functions-----------------------------------------------------
477// ------------------------------------------------------------------------------------------------------------
478
479template <class T>
481{
482 gsInfo<<"\t";
483 gsInfo<<std::setw(4)<<std::left<<"It.";
484 gsInfo<<std::setw(17)<<std::left<<"Res. F";
485 gsInfo<<std::setw(17)<<std::left<<"|dU|/|Du|";
486 gsInfo<<std::setw(17)<<std::left<<"dL/DL";
487 gsInfo<<std::setw(17)<<std::left<<"|U|";
488 gsInfo<<std::setw(17)<<std::left<<"L";
489 gsInfo<<std::setw(17)<<std::left<<"|DU|";
490 gsInfo<<std::setw(17)<<std::left<<"DL";
491 gsInfo<<std::setw(17)<<std::left<<"|dU|";
492 gsInfo<<std::setw(17)<<std::left<<"dL";
493 gsInfo<<std::setw(17)<<std::left<<"ds²";
494 gsInfo<<std::setw(17)<<std::left<<"|dU|²";
495 gsInfo<<std::setw(17)<<std::left<<"dL²";
496 gsInfo<<std::setw(17)<<std::left<<"Dmin";
497 gsInfo<<std::setw(17)<<std::left<<"note";
498 gsInfo<<"\n";
499
500 m_note = "";
501}
502
503template <class T>
505{
506 // if (!m_quasiNewton)
507 // {
508 computeStability(false);
509 // }
510 // else
511 // m_indicator = 0;
512
513 T A0 = math::pow(m_phi,2)*m_forcing.dot(m_forcing);
514
515 gsInfo<<"\t";
516 gsInfo<<std::setw(4)<<std::left<<m_numIterations;
517 gsInfo<<std::setw(17)<<std::left<<m_residueF;
518 gsInfo<<std::setw(17)<<std::left<<m_residueU;
519 gsInfo<<std::setw(17)<<std::left<<m_residueL;
520 gsInfo<<std::setw(17)<<std::left<<(m_U+m_DeltaU).norm();
521 gsInfo<<std::setw(17)<<std::left<<(m_L + m_DeltaL);
522 gsInfo<<std::setw(17)<<std::left<<m_DeltaU.norm();
523 gsInfo<<std::setw(17)<<std::left<<m_DeltaL;
524 gsInfo<<std::setw(17)<<std::left<<m_deltaU.norm();
525 gsInfo<<std::setw(17)<<std::left<<m_deltaL;
526 gsInfo<<std::setw(17)<<std::left<<this->distance(m_DeltaU,m_DeltaL);//math::pow(m_DeltaU.dot(m_DeltaU) + A0*math::pow(m_DeltaL,2.0),0.5);
527 gsInfo<<std::setw(17)<<std::left<<math::pow(m_DeltaU.norm(),2.0);
528 gsInfo<<std::setw(17)<<std::left<<A0*math::pow(m_DeltaL,2.0);
529 gsInfo<<std::setw(17)<<std::left<<m_indicator <<std::left << " (" <<std::left<< m_negatives<<std::left << ")";
530 gsInfo<<std::setw(17)<<std::left<<m_note;
531 gsInfo<<"\n";
532
533 m_note = "";
534}
535
536} // namespace gismo
Performs the Crisfield arc length method to solve a nonlinear equation system.
Definition gsALMCrisfield.h:30
void defaultOptions()
See gsALMBase.
Definition gsALMCrisfield.hpp:23
void predictor()
See gsALMBase.
Definition gsALMCrisfield.hpp:111
void computeLambdasModified()
Compute the load factors.
Definition gsALMCrisfield.hpp:251
void computeLambdaDET()
Compute the load factors.
Definition gsALMCrisfield.hpp:367
void initOutput()
See gsALMBase.
Definition gsALMCrisfield.hpp:480
void computeLambdaDOT()
Compute the load factors.
Definition gsALMCrisfield.hpp:402
void stepOutput()
See gsALMBase.
Definition gsALMCrisfield.hpp:504
void computeLambdasSimple()
Compute the load factors.
Definition gsALMCrisfield.hpp:213
void iteration()
See gsALMBase.
Definition gsALMCrisfield.hpp:74
void quasiNewtonPredictor()
See gsALMBase.
Definition gsALMCrisfield.hpp:56
void initiateStep()
See gsALMBase.
Definition gsALMCrisfield.hpp:103
void quasiNewtonIteration()
See gsALMBase.
Definition gsALMCrisfield.hpp:66
void initMethods()
See gsALMBase.
Definition gsALMCrisfield.hpp:41
void computeLambdaMU()
Compute the load factors.
Definition gsALMCrisfield.hpp:385
void iterationFinish()
See gsALMBase.
Definition gsALMCrisfield.hpp:194
void computeLambdasComplex()
Compute the load factors.
Definition gsALMCrisfield.hpp:308
void computeLambdas()
Compute the load factors.
Definition gsALMCrisfield.hpp:328
void getOptions()
See gsALMBase.
Definition gsALMCrisfield.hpp:31
void computeLambdasEta()
Compute the load factors.
Definition gsALMCrisfield.hpp:234
A vector with arbitrary coefficient type and fixed or dynamic size.
Definition gsVector.h:37
#define index_t
Definition gsConfig.h:32
#define gsInfo
Definition gsDebug.h:43
#define GISMO_ASSERT(cond, message)
Definition gsDebug.h:89
The G+Smo namespace, containing all definitions for the library.
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 <j> in the set ; by def...