1 //#************************************************************** 2 //# 3 //# filename: femath.h 4 //# 5 //# author: Gerstmayr Johannes 6 //# 7 //# generated: 15.05.97 8 //# description: Classes for linear and nonlinear algebra which is 9 //# thought to be used in finite element or similar applications 10 //# There are 2D and 3D and arbitrary size Vectors (Vector2D, Vector3D, Vector), 11 //# arbitrary size matrices (Matrix), and a nonlinear system (NumNLSys) 12 //# and a nonlinear solver (NumNLSolver) 13 //# remarks: Indizes run from 1 to n except in Vector3D/2D 14 //# 15 //# Copyright (c) 2003-2013 Johannes Gerstmayr, Linz Center of Mechatronics GmbH, Austrian 16 //# Center of Competence in Mechatronics GmbH, Institute of Technical Mechanics at the 17 //# Johannes Kepler Universitaet Linz, Austria. All rights reserved. 18 //# 19 //# This file is part of HotInt. 20 //# HotInt is free software: you can redistribute it and/or modify it under the terms of 21 //# the HOTINT license. See folder 'licenses' for more details. 22 //# 23 //# bug reports are welcome!!! 24 //# WWW: www.hotint.org 25 //# email: bug_reports@hotint.org or support@hotint.org 26 //#*************************************************************************************** 27 28 #ifndef MFEMATH__H 29 #define MFEMATH__H 30 31 #include <math.h> 32 #include "release_assert.h" 33 #include "solversettings_auto.h" 34 35 ///////////////////////// 36 //turns off critical asserts (faster): 37 //PG: Here we define the strategy for asserting femath code. 38 //We still have to discuss if this is the optimal setting. 39 //Meanwhile, I suggest: it's best to assert in debug-mode 40 //and if the preprocessor flag "__ASSERT_IN_RELEASE_MODE__" 41 //is set (e.g., in "MBSKernelLib\preprocessor_includes.h"). 42 //This is due to the fact, that some of the models may 43 //be executed in release-mode only. 44 #ifndef _DEBUG 45 #ifndef __ASSERT_IN_RELEASE_MODE__ 46 #define __QUICKMATH 47 #endif 48 #else 49 struct UserOutputInterface; 50 extern UserOutputInterface * global_uo; 51 #endif 52 ///////////////////////// 53 54 55 #include "mbs_interface.h" 56 57 #include "femathHelperFunctions.h" 58 59 typedef TArray <DVector*> DArray; 60 typedef TArray <IVector*> IArray; 61 typedef TArray <Vector*> VArray; 62 typedef TArray <Vector3D> AVector3D; 63 64 struct SuperMatrix; 65 //#include "..\SuperLU\supermatrix.h" 66 #include "..\SuperLU\SuperLUmain.h" //for external solver 67 68 #include "linalg.h" 69 #include "linalgeig.h" 70 #include "elementdata.h" //$AD 2011-03-24: new entry for use of parser in mathfunc 71 #include "parser.h" //$AD 2011-03-24: new entry for use of parser in mathfunc 72 #include "mathfunc.h" 73 #include <sstream> 74 75 #include "..\workingmodule\WorkingModuleBaseClass.h" 76 77 class NumNLSys; 78 79 80 class SparseJacMat 81 { 82 public: SparseJacMat()83 SparseJacMat(): J_vv_band_pivot(), J_vv(), J_zv(), J_vz(), J_zz() 84 , sparseinv(0) 85 { 86 LUcomputed = 0; 87 LUcomputed_bw = 0; 88 Init(); 89 }; 90 Init()91 void Init() 92 { 93 J_vv.Init(); 94 J_vv_band_pivot.Init(); 95 LUcomputed = 0; 96 LUcomputed_bw = 0; 97 useband = 0; 98 n_vv_written = 0; 99 100 sparseinv = 0; 101 102 } 103 ~SparseJacMat()104 ~SparseJacMat() 105 { 106 delete sparseinv; 107 } 108 109 void Apply(const Vector& R, Vector& d, NumNLSys* nls); 110 111 int Factorize(NumNLSys* nls); 112 113 void ApplyTransform(Vector& v, int mode, NumNLSys* nls); 114 115 SparseMatrix J_vv; 116 SparseMatrix J_zvS, J_vzS; 117 118 Matrix J_zv, J_vz, J_zz; 119 Matrix J_sch; 120 Matrix J_vv_band; 121 Matrix J_vv_band_aux; 122 IVector J_vv_band_pivot; 123 Matrix minv_test; 124 125 IVector LUindx; //temporary storage full LU decomposition 126 Vector LUvv; //temporary storage full LU decomposition 127 128 //++++++++++++++++++++++++++++ 129 // SUPERLU; PARDISO 130 131 SparseInverse *sparseinv; 132 //++++++++++++++++++++++++++++ 133 134 int lu; 135 int LUcomputed; 136 int LUcomputed_bw; 137 int n_vv_written; //set to 1 if state has been already written! 138 139 //for renumbering of band-diagonal part (rigid bodies!!!) 140 IVector resortvector; 141 int resortsize; 142 143 Vector temp; //temporary variables for Apply 144 Vector temp2; //temporary variables for Apply 145 Vector Rsort; //temporary variables for Apply 146 Vector vtrans; //temporary variable for ApplyTransform 147 148 Matrix mtemp1, mtemp2; //temporary variables for Factorize 149 Vector tempcol; //temporary variables for Factorize 150 151 int useband; 152 }; 153 154 class SaveJac 155 { 156 public: SaveJac()157 SaveJac(): oldjacmat(), oldsjacmat(), sparsejacmat() 158 { 159 oldjacmat.Init(); 160 oldsjacmat.Init(); 161 sparsejacmat.Init(); 162 Init(); 163 } 164 Init()165 void Init() 166 { 167 oldjac=0; //jacobimatrix existiert 168 oldjacsize=0; //gr��e 169 oldjacage=0; //alter 170 maxjacage=10000000;//10000000; //maximales alter 171 nlsinfo = 0; //nonlinear solve info (z.b. zeitschrittweite von timeint.cpp) 172 jacfailedcnt = 0; 173 lastbuilt = 0; 174 175 oldjacmat.SetSize(0,0); 176 oldsjacmat.SetSize(0,0); 177 sparsejacmat.Init(); 178 } 179 180 int oldjac; 181 int oldjacsize; 182 int oldjacage; 183 int maxjacage; 184 Matrix oldjacmat; 185 SparseMatrix oldsjacmat; 186 SparseJacMat sparsejacmat; 187 double nlsinfo; 188 int jacfailedcnt; 189 int lastbuilt; 190 }; 191 192 #include "options_class_auto.h" 193 194 class NumNLSolver; 195 196 //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 197 //Numerical nonlinear solver 198 //Needs a function which is solved for x f(x)=0 199 //Jacobion matrix is calculated by numerical differentiation 200 //$!YV 2012-11-29: inherits MBS interface, which is then implemented in the class hierarchy upwards 201 class NumNLSys : public WorkingModuleBaseClass, public MBS 202 { 203 204 public: NumNLSys()205 NumNLSys() {jaccol = 0; TIstages=1; jacfullnewton=0; fullnewtoncnt = 0; doresort = 0; hotint_options = new HOTINTOptions(this);}; ~NumNLSys()206 ~NumNLSys() { DeleteOptions(); } 207 virtual void NLF(const Vector& x, Vector& f)=0; //nonlinear function F(x) = 0 208 virtual void SetSolver(NumNLSolver* s); 209 virtual void Jacobian(Matrix& m, Vector& x); Jacobian(SparseMatrix & m,Vector & x)210 virtual void Jacobian(SparseMatrix& m, Vector& x) {}; 211 212 int NLS_GetJacCol() const; NLS_SetJacCol(int i)213 void NLS_SetJacCol(int i) {jaccol = i;} NLS_GetTIstages()214 int NLS_GetTIstages() const {return TIstages;} NLS_SetTIstages(int i)215 void NLS_SetTIstages(int i) {TIstages=i;}; NLS_IsJacFullNewton()216 int NLS_IsJacFullNewton() const {return jacfullnewton;} NLS_SetJacFullNewton(int i)217 void NLS_SetJacFullNewton(int i) {jacfullnewton=i;}; FullNewtonCnt()218 int FullNewtonCnt() const {return fullnewtoncnt;} FullNewtonCnt()219 int& FullNewtonCnt() {return fullnewtoncnt;} ResetNLSys()220 void ResetNLSys() {fullnewtoncnt = 0;} 221 GetResortVector()222 virtual const IVector& GetResortVector() const {return resortvector;} GetResortVector()223 virtual IVector& GetResortVector() {return resortvector;} GetResortSize()224 virtual int GetResortSize() const {return resortsize;}; SetResortSize(int s)225 virtual void SetResortSize(int s) {resortsize = s;}; 226 227 virtual int UseSparseSolver() const; 228 virtual int SolveUndeterminedSystem() const; 229 virtual int& SolveUndeterminedSystem(); 230 virtual double EstimatedConditionNumber() const; 231 virtual double& EstimatedConditionNumber(); 232 virtual int UseSparseJac() const; 233 virtual int UseSuperLU() const; 234 TransformJacApply()235 virtual int TransformJacApply() const {return 0;} ApplyTransform(const Vector & v,Vector & Av,int mode)236 virtual void ApplyTransform(const Vector& v, Vector& Av, int mode) {}; //compute Av=A^T*v in mode==0 and Av=A*v in mode==1 237 238 //[I|D|T]Options 239 virtual void SetIOption(int index, int data); 240 virtual const int& GetIOption(int index) const; 241 virtual int& GetIOption(int index); 242 virtual void SetDOption(int index, double data); 243 virtual const double& GetDOption(int index) const; 244 virtual double& GetDOption(int index); 245 246 virtual void SetTOption(int index, const char* data); 247 virtual const char* GetTOption(int index) const; 248 virtual void InitializeOptions(); 249 //+++++++++++++++++++++++++++++++++++++++++++++++++ 250 // implemented in auto - generated file "options_auto.h" 251 virtual void InitializeOptionsAuto(); 252 //+++++++++++++++++++++++++++++++++++++++++++++++++ 253 virtual void DeleteOptions(); 254 GetOptions()255 virtual HOTINTOptions* GetOptions() { return hotint_options; } GetOptions()256 virtual const HOTINTOptions* GetOptions() const { return hotint_options; } 257 258 ///////////////////////////// 259 //#output details of computation: (changes in the 'Edit All Options' window take effect immediately, since the solset.* values are returned here!) 260 261 // definition of user output level for detailed solver logs UO_LVL_SolverDetails()262 UO_MSGLVL UO_LVL_SolverDetails() const { return UO_LVL_all; } 263 264 // solver details are only printed to logfile, if file_output_level is set greater or equal to UO_LVL_SolverDetails() 265 // solver details are only printed to window, if output_level is set greater or equal to UO_LVL_SolverDetails() SolverPrintsDetailedOutput()266 bool SolverPrintsDetailedOutput() 267 { 268 return 269 ( 270 SolverUO().GetGlobalFileMessageLevel() >= UO_LVL_SolverDetails() || 271 SolverUO().GetGlobalMessageLevel() >= UO_LVL_SolverDetails() 272 ) 273 && 274 AnySolverOutputChecked(); 275 } 276 AnySolverOutputChecked()277 bool AnySolverOutputChecked() 278 { 279 return 280 GetOptions()->LoggingOptions()->SolverGeneralInformation() || 281 GetOptions()->LoggingOptions()->SolverNewtonIterationJacobiMatrix() || 282 GetOptions()->LoggingOptions()->SolverNewtonIterationResidualVector() || 283 GetOptions()->LoggingOptions()->SolverNewtonIterationSolutionVector() || 284 GetOptions()->LoggingOptions()->SolverPostNewtonIterationDataVector() || 285 GetOptions()->LoggingOptions()->SolverStepSolutionVector() || 286 GetOptions()->LoggingOptions()->SolverStepSolutionVectorIncrement(); 287 } 288 SolverUO()289 UserOutput& SolverUO() { return UOfull(UO_LVL_SolverDetails()); } SolverUO()290 UserOutput& SolverUO() const { return UOfull(UO_LVL_SolverDetails()); } 291 292 //virtual UserOutput& UO(int message_level = UO_LVL_all) {uo.SetLocalMessageLevel(message_level); return uo;}; //(AD) 293 virtual UserOutput& UOfull(int message_level = UO_LVL_all, int output_prec = -1) const; //$ AD 2011-02 output_prec 294 virtual UserOutputInterface& UO(int message_level = UO_LVL_all, int output_prec = -1) const { return UOfull(message_level,output_prec); } 295 296 //option members: 297 int ioptions[301]; 298 double doptions[301]; 299 char* toptions[301]; 300 301 302 private: 303 int jaccol; 304 int jacfullnewton; 305 int TIstages; 306 NumNLSolver* solver; 307 int fullnewtoncnt; 308 309 Vector f1; //temporary variables for Jacobian 310 Vector f2; //temporary variables for Jacobian 311 312 protected: 313 int resortsize; 314 int doresort; 315 IVector resortvector; 316 HOTINTOptions* hotint_options; 317 }; 318 319 320 class NumNLSolver : public NumSolverInterface 321 { 322 SolverSettings* solset; //this part of the solversettings are referenced to the settings, which can be online modified by user; do only use solset for options which should show 323 //immediate effect of change, e.g. the logging of newton iterations or errors; DO NOT USE E.G. FOR NEWTON PARAMETERS (e.g. TOLERANCES)==>THIS LEADS TO ARBITRARY RESULTS 324 public: NumNLSolver()325 NumNLSolver():iv() {}; NumNLSolver(NumNLSys * nlsi,SolverSettings * solsetI)326 NumNLSolver(NumNLSys* nlsi, SolverSettings* solsetI):iv() 327 { 328 nls = nlsi; 329 solset = solsetI; 330 nls->SetSolver(this); 331 332 relativeaccuracy = 1e-8; 333 maxmodnewtonsteps=50; 334 maxrestartnewtonsteps=25; 335 maxfullnewtonsteps=20; 336 numdiffepsi = 1e-8; 337 newtonits = 0; 338 trustregion = 0; 339 //trustregionitmax = 6; 340 trustregiondiv = 1./10.; 341 output=0; 342 modifiednewton = 1; // Sets solving-method to modified Newton 343 symmetricjacobian = 1; 344 nonlinsolveinfo = 0; 345 stopmnr = 0; 346 usesparsesolver = 0; 347 solveundeterminedsystem = 0; 348 estimatedconditionnumber = 1e12; 349 350 contractivity = 0; 351 352 jaccount=0; 353 fulljaccnt = 0; 354 355 jac_condnum=-1; 356 error_msg = ""; 357 358 bandsize = 0; 359 360 //DestroyOldJacs(); 361 //sjac[0].oldjacmat = Matrix(); 362 //sjac[1].oldjacmat = Matrix(); 363 364 //log.SetAllInactive(); 365 } 366 ResetSolver()367 void ResetSolver() 368 { 369 relativeaccuracy = 1e-8; 370 maxmodnewtonsteps=50; 371 maxrestartnewtonsteps=25; 372 maxfullnewtonsteps=20; 373 numdiffepsi = 1e-8; 374 newtonits = 0; 375 trustregion = 0; 376 //trustregionitmax = 6; 377 trustregiondiv = 1./10.; 378 output=0; 379 modifiednewton = 1; // Sets solving-method to modified Newton 380 symmetricjacobian = 1; 381 nonlinsolveinfo = 0; 382 stopmnr = 0; 383 384 contractivity = 0; 385 386 jaccount=0; 387 fulljaccnt = 0; 388 389 jac_condnum=-1; 390 error_msg = ""; 391 392 bandsize = 0; 393 394 DestroyOldJacs(); 395 //sjac[0].oldjacmat = Matrix(); 396 //sjac[1].oldjacmat = Matrix(); 397 sjac[0].Init(); 398 sjac[1].Init(); 399 400 } 401 402 //Newton Solver 403 int NLSolve(Vector& x0); 404 405 int Factorize(Matrix& minv, SparseMatrix& sminv, SparseJacMat& sparseminv); 406 407 //set res=Jac^(-1)*f if jac exists, else return 0 408 int ApplyJac(const Vector& f, Vector& res); 409 ChooseJac()410 int ChooseJac() 411 { 412 if (sjac[0].oldjac && RelApproxi(sjac[0].nlsinfo,nonlinsolveinfo,2e-2)) 413 { 414 return 1; 415 } 416 else if (sjac[1].oldjac && RelApproxi(sjac[1].nlsinfo,nonlinsolveinfo,2e-2)) 417 { 418 return 2; 419 } 420 else return 0; 421 } 422 423 //virtual void Jacobian(Matrix& m, Vector& x); 424 GetNewtonIts()425 int GetNewtonIts() const {return newtonits;}; 426 ModifiedNewton()427 int ModifiedNewton() const {return modifiednewton;}; ModifiedNewton()428 int& ModifiedNewton() {return modifiednewton;}; 429 AbsoluteAccuracy()430 double AbsoluteAccuracy() const {return absoluteaccuracy;}; AbsoluteAccuracy()431 double& AbsoluteAccuracy() {return absoluteaccuracy;}; 432 RelativeAccuracy()433 double RelativeAccuracy() const {return relativeaccuracy;}; RelativeAccuracy()434 double& RelativeAccuracy() {return relativeaccuracy;}; 435 436 //int& MaxNewtonSteps() {return maxnewtonsteps;}; //$ PG 2013-9-19: who changed that??? 437 //int MaxNewtonSteps() const {return MaxModNewtonSteps()+MaxRestartNewtonSteps()+MaxFullNewtonSteps();}; //$ PG 2013-9-19: does not make any sense in case of full newton! MaxNewtonSteps()438 int MaxNewtonSteps() const 439 { 440 if (!ModifiedNewton()) 441 { 442 return MaxFullNewtonSteps(); 443 } 444 445 return MaxModNewtonSteps()+MaxRestartNewtonSteps()+MaxFullNewtonSteps(); 446 } 447 MaxModNewtonSteps()448 int MaxModNewtonSteps() const {return maxmodnewtonsteps;}; MaxModNewtonSteps()449 int& MaxModNewtonSteps() {return maxmodnewtonsteps;}; 450 MaxRestartNewtonSteps()451 int MaxRestartNewtonSteps() const {return maxrestartnewtonsteps;}; MaxRestartNewtonSteps()452 int& MaxRestartNewtonSteps() {return maxrestartnewtonsteps;}; 453 MaxFullNewtonSteps()454 int MaxFullNewtonSteps() const {return maxfullnewtonsteps;}; MaxFullNewtonSteps()455 int& MaxFullNewtonSteps() {return maxfullnewtonsteps;}; 456 NumDiffepsi()457 double NumDiffepsi() const {return numdiffepsi;}; NumDiffepsi()458 double& NumDiffepsi() {return numdiffepsi;}; 459 TrustRegion()460 int TrustRegion() const {return trustregion;}; TrustRegion()461 int& TrustRegion() {return trustregion;}; 462 463 //int TrustRegionItMax() const {return trustregionitmax;}; 464 //int& TrustRegionItMax() {return trustregionitmax;}; 465 TrustRegionDiv()466 double TrustRegionDiv() const {return trustregiondiv;}; TrustRegionDiv()467 double& TrustRegionDiv() {return trustregiondiv;}; 468 Output()469 int Output() const {return output;}; Output()470 int& Output() {return output;}; 471 SymmetricJacobian()472 int SymmetricJacobian() const {return symmetricjacobian;}; SymmetricJacobian()473 int& SymmetricJacobian() {return symmetricjacobian;}; 474 NLSolveInfo()475 double NLSolveInfo() const {return nonlinsolveinfo;}; NLSolveInfo()476 double& NLSolveInfo() {return nonlinsolveinfo;}; 477 Contractivity()478 double Contractivity() const {return contractivity;} 479 GetJac()480 const Matrix& GetJac() const 481 { 482 if (sjac[0].lastbuilt) return sjac[0].oldjacmat; 483 else return sjac[1].oldjacmat; 484 } DestroyOldJac(double info)485 void DestroyOldJac(double info) 486 { 487 if (info == sjac[0].nlsinfo) 488 { 489 sjac[0].oldjac = 0; sjac[0].nlsinfo = -1; sjac[0].lastbuilt = 0; sjac[1].lastbuilt = 1; 490 } 491 else if (info == sjac[1].nlsinfo) 492 { 493 sjac[1].oldjac = 0; sjac[1].nlsinfo = -1; sjac[0].lastbuilt = 1; sjac[1].lastbuilt = 0; 494 } 495 496 } DestroyOldJacs()497 void DestroyOldJacs() 498 { 499 sjac[0].oldjac = 0; sjac[0].nlsinfo = -1; sjac[0].lastbuilt = 0; 500 sjac[1].oldjac = 0; sjac[1].nlsinfo = -1; sjac[1].lastbuilt = 0; 501 } 502 ModifiedNewtonActive()503 int ModifiedNewtonActive() const {return ModifiedNewton() && (!stopmnr);} 504 GetJacCount()505 int GetJacCount() const {return jaccount;}; GetFullJacCnt()506 int GetFullJacCnt() const {return fulljaccnt;}; GetJacCondnum()507 double GetJacCondnum() const {return jac_condnum;}; GetErrorMsg()508 const mystr& GetErrorMsg() const {return error_msg;}; GetErrorMsg()509 mystr& GetErrorMsg() {return error_msg;}; 510 GetBandSize()511 int GetBandSize() const {return bandsize;} SetBandSize(int bs)512 void SetBandSize(int bs) {bandsize = bs;} 513 UseSparseSolver()514 int UseSparseSolver() const { return usesparsesolver; } UseSparseSolver()515 int& UseSparseSolver() { return usesparsesolver; } SolveUndeterminedSystem()516 int SolveUndeterminedSystem() const { return solveundeterminedsystem; } SolveUndeterminedSystem()517 int& SolveUndeterminedSystem() { return solveundeterminedsystem; } EstimatedConditionNumber()518 double EstimatedConditionNumber() const { return estimatedconditionnumber; } EstimatedConditionNumber()519 double& EstimatedConditionNumber() { return estimatedconditionnumber; } 520 521 522 //print detailed output concerning NLSolve computations? (time critical!) SolverPrintsDetailedOutput()523 bool SolverPrintsDetailedOutput() const { return nls->SolverPrintsDetailedOutput(); } SolverUO()524 UserOutput& SolverUO() { return nls->SolverUO(); } SolverUO()525 UserOutput& SolverUO() const { return nls->SolverUO(); } GetOptions()526 HOTINTOptions* GetOptions() { return nls->GetOptions(); } 527 528 529 private: 530 531 //$ PG 2013-11-6: assemble Jacobian (and print some information), set jaccount++ 532 void AssembleJacobian(Matrix* minv, SparseMatrix* sminv, Vector& x0); 533 534 //$ PG 2013-10-17: calculate condition number of jacobi matrix, IS SLOW, should only be used for designing elements, constraints, for adjusting models, or for general debugging purpose 535 void OutputConditionNumber(const SparseMatrix& M) const; 536 void OutputConditionNumber(const Matrix& M) const; 537 void OutputJacobian(const SparseMatrix& M) const; 538 void OutputJacobian(const Matrix& M) const; 539 540 int newtonits; 541 NumNLSys* nls; 542 543 SaveJac sjac[2]; 544 545 //solver settings 546 int modifiednewton; //indicates modified Newton algorithm is used 547 double absoluteaccuracy; 548 double relativeaccuracy; 549 double numdiffepsi; 550 int maxmodnewtonsteps, maxrestartnewtonsteps, maxfullnewtonsteps; 551 int stopmnr; 552 553 int trustregion; //turn on or off 554 //int trustregionitmax; 555 double trustregiondiv; 556 557 int output; 558 559 int symmetricjacobian; 560 561 double nonlinsolveinfo; 562 int usesparsesolver; 563 int solveundeterminedsystem; 564 double estimatedconditionnumber; 565 566 int jaccount; 567 double contractivity; 568 int fulljaccnt; 569 double jac_condnum; 570 mystr error_msg; 571 572 //sparse: 573 int bandsize; //size of leftupper system (first n unknowns) which has small bandwidth, esp. for MBS 574 575 576 Vector x0start, f, xd, hv; //temporary variables in NLsolve 577 IVector iv; //temporary variables in NLsolve 578 579 }; 580 581 582 583 #endif 584