1 //#************************************************************** 2 //# 3 //# filename: timeint.h 4 //# 5 //# author: Gerstmayr Johannes 6 //# 7 //# generated: 09.06.99 8 //# description: Class for implicit and explicit Runge Kutta Time integration, 9 //# variable stepsize and arbitrary order 10 //# remarks: The file tableaus.txt must be provided in Project/Release or Debug !!! 11 //# 12 //# Copyright (c) 2003-2013 Johannes Gerstmayr, Linz Center of Mechatronics GmbH, Austrian 13 //# Center of Competence in Mechatronics GmbH, Institute of Technical Mechanics at the 14 //# Johannes Kepler Universitaet Linz, Austria. All rights reserved. 15 //# 16 //# This file is part of HotInt. 17 //# HotInt is free software: you can redistribute it and/or modify it under the terms of 18 //# the HOTINT license. See folder 'licenses' for more details. 19 //# 20 //# bug reports are welcome!!! 21 //# WWW: www.hotint.org 22 //# email: bug_reports@hotint.org or support@hotint.org 23 //#*************************************************************************************** 24 25 26 #ifndef TIMEINT__H 27 #define TIMEINT__H 28 29 #include "mystring.h" 30 #include "timeintlog.h" 31 #include "femath.h" 32 33 //# maximum stages for static variables in K-Version timeint 34 #define global_maxstages 21 35 36 37 38 class IRK_Tableau 39 { 40 public: IRK_Tableau()41 IRK_Tableau() {}; 42 43 mystr name; 44 mystr info; 45 int nstages; 46 int ODE_order; 47 int DAE_order_y[10]; //DAE_order[1]=index1 order, DAE_order[2]=index2 order, DAE_order[3]=index3 order, ... 48 int DAE_order_z[10]; //DAE_order[1]=index1 order, DAE_order[2]=index2 order, DAE_order[3]=index3 order, ... 49 50 int Ainvertable; 51 int implicit; 52 53 Vector b; 54 Vector c; 55 Matrix A; 56 Matrix A2; 57 Matrix Ainv; 58 Vector bAinv; ///b^T * A^(-1) 59 }; 60 61 62 63 64 65 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 66 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 67 class TimeIntVars 68 { 69 public: TimeIntVars()70 TimeIntVars() {Init();} Init()71 virtual void Init() 72 { 73 start_clock_time = 0; 74 timetogo = 0; 75 lastdraw = 0; 76 lastprintres = 0; 77 last_showstatustext = 0; 78 timetogo2 = 0; 79 outputlevel = 3; //2 80 maxnewtonit=0; 81 } 82 83 //time, drawing, log messages: 84 double start_clock_time; //stores when clock was started 85 double timetogo; //estimated time to go 86 double lastdraw; //time when last was drawn 87 double lastprintres; //time when last was printed results 88 double last_showstatustext; //time when last was printed results 89 double timetogo2; //damping term for time to go estimation 90 int outputlevel; //output level, how much is printed per step 91 92 //temporary solver variables: 93 int maxnewtonit; //maximum number of newton iterations since last printing 94 }; 95 96 97 98 //structure of data and state in one time step 99 struct TIstepsave 100 { 101 Vector x0; //state variables 102 Vector data; //data variables 103 Vector k0; //state variables in differentiated form 104 double time; 105 }; 106 107 108 class TimeInt : public NumNLSys 109 { 110 public: 111 TimeInt(); 112 ~TimeInt()113 virtual ~TimeInt() 114 { 115 for (int i=1; i <= tableaus.Length(); i++) 116 { 117 delete tableaus.Elem(i); 118 } 119 } 120 //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 121 //#main functions for time integration 122 //#x = [u v q z] // \dot u = v, M (TIx0) \dot v = F2(TIx0,t), \dot q = F2(TIx0,t), G(TIx0) = 0 123 //#len(u) = len(v), [len(u)=0 | len(q) = 0], len(z) arbitrary PrecomputeEvalFunctions()124 virtual void PrecomputeEvalFunctions() {}; //$ PG 2012-1-15: precomputation for EvalF(), EvalF2(), EvalG(), EvalM(), and also for CqTLambda() in case of constraints EvalF(const Vector & x,Vector & f,double t)125 virtual void EvalF(const Vector& x, Vector& f, double t) {}; //*first order equations: \dot q=f(q,t), return vector: len(q) EvalM(const Vector & x,Matrix & m,double t)126 virtual void EvalM(const Vector& x, Matrix& m, double t) {}; //*evaluate mass matrix M(TIx0), return matrix: len(u)*len(u) EvalM(const Vector & x,SparseMatrix & m,double t)127 virtual void EvalM(const Vector& x, SparseMatrix& m, double t) {}; //*evaluate mass matrix M(TIx0), return matrix: len(u)*len(u) SetBWM(int i)128 virtual void SetBWM(int i) {bandwidthm = i;} EvalF2(const Vector & x,Vector & f,double t)129 virtual void EvalF2(const Vector& x, Vector& f, double t) {}; //*second order equations: M \ddot u = F2(u,\dot u,t), len(u) EvalMinvF2(const Vector & x,Vector & f,double t)130 virtual void EvalMinvF2(const Vector& x, Vector& f, double t) {}; //*second order equations: M \ddot u = F2(u,\dot u,t), len(u) 131 //virtual void EvalF2(const Vector& x, Vector& f, double t, int& minind, int& maxind) {}; //only for jacobian, provide also min and max entry EvalF2(const Vector & x,SparseVector & f,double t,IVector & rowind,IVector & clearind,IVector & elems)132 virtual void EvalF2(const Vector& x, SparseVector& f, double t, 133 IVector& rowind, IVector& clearind, IVector& elems) {}; //only for jacobian, rowind and clearind for speedup, rowind must be initialized with zeros and l=f.Length() EvalFelastic(const Vector & x,Vector & f,double t)134 virtual void EvalFelastic(const Vector& x, Vector& f, double t) {}; 135 EvalG(const Vector & x,Vector & f,double t)136 virtual void EvalG(const Vector& x, Vector& f, double t) {}; //*evaluate constraints: len(z) EvalG(const Vector & x,SparseVector & f,double t,IVector & rowind,IVector & clearind)137 virtual void EvalG(const Vector& x, SparseVector& f, double t, IVector& rowind, IVector& clearind) {}; //only for jacobian, rowind and clearind for speedup, rowind must be initialized with zeros and l=f.Length() 138 GetSecondOrderSize()139 virtual int GetSecondOrderSize() const {return 0;}; //*size of second order equations, len(u) GetSecondOrderSize_RS()140 virtual int GetSecondOrderSize_RS() const {return GetSecondOrderSize();}; //*size of second order equations, len(u), after resorting!!! GetFirstOrderSize()141 virtual int GetFirstOrderSize() const {return 0;}; //*size of first order equations GetImplicitSize()142 virtual int GetImplicitSize() const {return 0;}; //* len(z) GetSystemSize()143 virtual int GetSystemSize() const {return GetSecondOrderSize()*2 144 +GetFirstOrderSize()+GetImplicitSize();}; //*len(z)+len(u)+len(v)+len(q) GetDataSize()145 virtual int GetDataSize() const {return 0;} //size of data variables for each time step of system 146 147 148 virtual void NLF(const Vector& x, Vector& f); //* 149 virtual int NonlinSubStep(); //* 150 virtual void FinishStep(); //* 151 virtual void TIDrawAndStore(); //* GetEnergy()152 virtual double GetEnergy() const {return 0;}; Initialize()153 virtual void Initialize() {}; //initialize MBS system SetInitialConditions()154 virtual void SetInitialConditions() {}; //initial conditions and some initializations of MBS 155 virtual void InitFirst(); //initialize at program start when UserInterface is linked, before config file InitializeAfterConfigLoaded()156 virtual void InitializeAfterConfigLoaded() {}; 157 //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 158 //#internal functions for simulation 159 void LoadTableaus(const mystr& tableauname); // nur filename mit tableaus 160 int GetTableau(mystr tabname, int stages); //tableau-nummer ausw�hlen 161 int LoadTableauList(mystr tabname, 162 TArray<int>& list, 163 int minstage, int maxstage); GetActTableau()164 IRK_Tableau* GetActTableau() const {return tableaus[act_tableau];} 165 166 virtual void TIInit(); //* 167 virtual int FullStep(); //* 168 virtual int FirstStep(); 169 virtual int GeneralIStep(); 170 GetCharacteristicSol()171 virtual double GetCharacteristicSol() const {return TIx0(1);} 172 173 virtual void PrintTimingList(); 174 virtual void ResetComputation(); 175 virtual int Integrate();//JG Integrate(const SolverSettings& solversettings); 176 177 virtual int IntegrateStep(); //integrate one step (adaptive or constant stepsize) 178 virtual int StaticComputation(); //try to solve static problem 179 180 //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 181 virtual int ExplicitIntegration(); //do explicit integration, return 1 if success 182 virtual int AdaptiveExplicitStep(); //perform one adaptive explicit step, return 1 if success 183 virtual int ExplicitStep(); //do one explicit nonlinear step, return 1 if success 184 virtual int GeneralEStep(); //Do one general explicit step 185 //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 186 187 //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 188 //# Computation Steps with changing SolverOptions... 189 //Time mapping functions ParseStepEndTimes()190 virtual int ParseStepEndTimes() {return 0;}; // parses all EDC for ComputationSteps entries - implemented in MBS DoubleCheckWithLoadSettings()191 virtual int DoubleCheckWithLoadSettings() {return 0;}; 192 193 TArray<double> CSEndTimes; // array the holds all endtimes of the computation steps GetTMaxCompSteps()194 virtual double GetTMaxCompSteps() {return CSEndTimes.Last();}; GetEndTimeCompStep(int stepnumber)195 virtual double GetEndTimeCompStep(int stepnumber) 196 { 197 if (stepnumber <= CSEndTimes.Length()) return CSEndTimes(stepnumber); 198 else return solset.endtime; // in case last computation step finishes BEFORE end of simulation 199 }; 200 201 int computationstepnumber; // number of the computation step ("current") ComputationStepNumber()202 int& ComputationStepNumber() {return computationstepnumber;}; // quick access to the computationstepnumber (variable, no recalculaiton) 203 double computationsteptime; // relative time within the current computation step ("percentage") ComputationStepTime()204 double& ComputationStepTime() {return computationsteptime;}; // quick access to the computaionsteptime (variable, no recalculation) 205 206 virtual int GetCSNumber(double globaltime); // calculates number of the Computation Step for given global time // returns -1 it globaltime < 0 or globaltime > last step_end_time 207 virtual double GetCSTime(double globaltime); // calculates relative time in corresponding step rv(0..1] // returns 0. only for globaltime == 0 or globaltime > last step_end_time 208 virtual int IsComputationStepFinished(double globaltime); // returns 1 if globaltime is any step_end_time 209 210 //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 211 //# Computation Steps with changing SolverOptions... 212 //replace solver options ApplyComputationStepSettings(int stepnumber)213 virtual int ApplyComputationStepSettings(int stepnumber) {return 0;} // changes the solver options - implemented in MBS ComputationStepFinished()214 virtual void ComputationStepFinished() {}; //function is called when a computation step is finished 215 216 //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 217 virtual void StepPrintAndDraw(int rv); //print step information and draw results each step 218 SetComputationSolverOptions()219 virtual void SetComputationSolverOptions() {}; GetTMaxLoadStepsLength()220 virtual int GetTMaxLoadStepsLength() const {return 0;} //overloaded in mbs.h 221 //Compute weights for Lagrangian interpolation, tau=0..1 222 double LagrangeWeight(int l, double tau, const Vector& c) const; 223 224 virtual int FullAdaptiveStep(); 225 StartTimeStep()226 virtual void StartTimeStep() {}; //function is called when computation of time step is started EndTimeStep()227 virtual void EndTimeStep() {}; //function is called when computation of time step is started 228 PostNewtonStep(double)229 virtual double PostNewtonStep(double ) {return 0;}; PostprocessingStep()230 virtual void PostprocessingStep() {}; 231 232 //#save and restore inner variables for calculating the error, restarting, ... SaveState()233 virtual void SaveState() {}; RestoreState()234 virtual void RestoreState() {}; 235 virtual void TISaveState(TIstepsave& stepsave); 236 virtual void TIRestoreState(const TIstepsave& stepsave); 237 virtual void LocalJacobianF(Matrix& m, Vector& x); //Compute (df/dx) 238 virtual void LocalJacobianF2(Matrix& m, Vector& x); //Compute (df2/dx) 239 virtual void LocalJacobianF2(SparseMatrix& m, Vector& x); //Compute (df2/dx) 240 virtual void LocalJacobianG(Matrix& m, Vector& x); //Compute (dg/dx) 241 virtual void LocalJacobianG(SparseMatrix& m, Vector& x); //Compute (df2/dx) 242 virtual void LocalJacobianM(Matrix& m, Vector& x); //Compute (d(M*TIkv)/dx) 243 virtual void Jacobian(Matrix& m, Vector& x); //Compute approximate Jacobian (if modified newton) or call base function 244 virtual void Jacobian(SparseMatrix& m, Vector& x); //Compute approximate Jacobian (if modified newton) or call base function 245 virtual void StaticJacobian(SparseMatrix& m, Vector& x); //Compute approximate Jacobian for static case (if modified newton) or call base function 246 virtual void FullJacobian(Matrix& m, Vector& x); //Compute approximate Jacobian (if modified newton) or call base function ForceJacobianRecomputation()247 virtual void ForceJacobianRecomputation() {numsol.DestroyOldJacs();} //Destroy actual jacobians and compute new jacobian, e.g. if stiffness or configuration of system changes has been changed significantly 248 249 250 virtual void ComputeStiffnessAndDampingMatrix(SparseMatrix& k, SparseMatrix& d, Vector& x); //Compute (df2/dx),(df2/dxp), directly fill into final stiffness and damping matrix 251 virtual void ComputeGyroscopicMatrix(SparseMatrix& gy); // $ MSax 2013-07-25 : added 252 StaticJacobianF2(SparseMatrix & m,Vector & x)253 virtual void StaticJacobianF2(SparseMatrix& m, Vector& x) {}; //Compute (df2/dx), directly fill into final stiffness matrix 254 ComputeSparseMatrixUsage(IVector & usage_per_dof)255 virtual void ComputeSparseMatrixUsage(IVector& usage_per_dof) {}; //compute entries of each row in sparse matrix 256 UseSparseMass()257 virtual int UseSparseMass() const {return UseSparseJac();} MaxSparseBandwidth()258 virtual int MaxSparseBandwidth() const {return 32;} //default value SetReducedBandsize(int reducedbandsizeI)259 virtual void SetReducedBandsize(int reducedbandsizeI) {reducedbandsize = reducedbandsizeI;} DoStaticComputation()260 virtual int DoStaticComputation() const {return solset.dostaticcomputation;} DoStaticComputation()261 virtual int& DoStaticComputation() {return solset.dostaticcomputation;} DoImplicitIntegration()262 virtual int DoImplicitIntegration() const {return solset.doimplicitintegration;} DoImplicitIntegration()263 virtual int& DoImplicitIntegration() {return solset.doimplicitintegration;} LoadFact()264 virtual double LoadFact() {return loadfact;} 265 266 267 //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 268 //#Access to solution and computation 269 270 //#starting the simulation: 271 virtual void SetStartVector(const Vector& x0); 272 //#get characteristic value for error; this error is attempted to be lower than the desired accuracy GetError()273 virtual double GetError() {return TIx0.GetNorm();}; 274 TIFinished()275 virtual void TIFinished() {TIfinished = 1;} TIFinishedWithError()276 virtual void TIFinishedWithError() {TIfinished = -1;} GetTIFinished()277 virtual const int& GetTIFinished() const {return TIfinished;} GetTIFinished()278 virtual int& GetTIFinished() {return TIfinished;} GetTIit()279 virtual int GetTIit() const {return TIit;} 280 GetNewtonItSum()281 virtual int GetNewtonItSum() const {return log.TInewtonitsum;} GetNewtonIt()282 virtual int GetNewtonIt() const {return log.TInewtonit;} SetNewtonItSum(int ii)283 virtual void SetNewtonItSum(int ii) {log.TInewtonitsum = ii;} SetNewtonIt(int ii)284 virtual void SetNewtonIt(int ii) {log.TInewtonit = ii;} 285 //virtual Vector GetNonlinSol() {return Vector(0);} //erase GetNonlinIts()286 virtual int GetNonlinIts() const {return log.TInonlinit;} 287 GetSolVector()288 virtual const Vector& GetSolVector() const {return TIx0;} //return actual states GetSolVector()289 virtual Vector& GetSolVector() {return TIx0;} SetSolVector(const Vector & sol)290 virtual void SetSolVector(const Vector& sol) {TIx0 = sol;} GetSol(int i)291 virtual const double& GetSol(int i) const {return TIx0(i);} 292 GetDataVector()293 virtual const Vector& GetDataVector() const {return TIdata;} //return non-state data vector (plastic strains, contact conditions, discrete states, etc.) GetDataVector()294 virtual Vector& GetDataVector() {return TIdata;} SetDataVector(const Vector & data)295 virtual void SetDataVector(const Vector& data) {TIdata = data; TIdrawdata = data;} //set non-state data vector (plastic strains, contact conditions, discrete states, etc.) 296 297 //$EK 2012-09-17: return non-state data draw vector (plastic strains, contact conditions, discrete states, etc.) GetDataDrawVector()298 virtual const Vector& GetDataDrawVector() const {return TIdrawdata;} GetDataDrawVector()299 virtual Vector& GetDataDrawVector() {return TIdrawdata;} 300 GetDataAct(int i)301 virtual const double& GetDataAct(int i) const {return TIdata(i);}; GetDataAct(int i)302 virtual double& GetDataAct(int i) {return TIdata(i);}; 303 GetDataDraw(int i)304 virtual const double& GetDataDraw(int i) const {return TIdrawdata(i);}; GetDataDraw(int i)305 virtual double& GetDataDraw(int i) {return TIdrawdata(i);}; 306 GetLastSolVector()307 virtual const Vector& GetLastSolVector() const {return TIlaststep_state.x0;} //return solution after end of last step == beginning of this step GetLastNLItSolVector()308 virtual const Vector& GetLastNLItSolVector() const {return TIlastnonlinit_state.x0;} //return solution after last nonlinear iteration 309 GetLastDataVector()310 virtual Vector& GetLastDataVector() {return TIlaststep_state.data;} //return data vector after end of last step == beginning of this step GetLastDataVector()311 virtual const Vector& GetLastDataVector() const {return TIlaststep_state.data;} //return data vector after end of last step == beginning of this step GetTempDataVector()312 virtual Vector& GetTempDataVector() {return TItemp_state.data;} //return data vector after StartTimeStep() GetTempDataVector()313 virtual const Vector& GetTempDataVector() const {return TItemp_state.data;} //return data vector after StartTimeStep() GetLastNLItDataVector()314 virtual Vector& GetLastNLItDataVector() {return TIlastnonlinit_state.data;} //return data vector after last nonlinear iteration GetLastNLItDataVector()315 virtual const Vector& GetLastNLItDataVector() const {return TIlastnonlinit_state.data;} //return data vector after last nonlinear iteration GetNumberOfSolutionDataSteps()316 virtual const int& GetNumberOfSolutionDataSteps() const { return TInumber_of_solution_data_steps; } //return number of stored soltuion data steps (which is particularly used for communication with data manager of WCDriver) SetNumberOfSolutionDataSteps(const int n_steps)317 virtual void SetNumberOfSolutionDataSteps(const int n_steps) { TInumber_of_solution_data_steps = n_steps; } //set number of stored soltuion data steps (which is particularly used for communication with data manager of WCDriver) 318 GetVelocityAndAccelerationVector()319 virtual const Vector& GetVelocityAndAccelerationVector() const {return TIk0;} //for 2nd order components, 1..ss // $ MSax 2013-07-16 : renamed from GetAccelerationVector to GetVelocityAndAccelerationVector GetVelocityAndAcceleration(int i)320 virtual const double& GetVelocityAndAcceleration(int i) const {return TIk0(i);} //1..ss: 1..sos:velocities, sos+1 .. 2*sos:accelerations, 2*sos+1 .. 2*sos+es:first order velocities, 2*sos+es+1..:algebraic variables // $ MSax 2013-07-16 : renamed from GetAcceleration to GetVelocityAndAcceleration 321 GetStepSize()322 virtual double GetStepSize() const {return TIstep;} SetStepSize(double s)323 virtual void SetStepSize(double s) {TIstepnew = s;} GetStepSizeNew()324 virtual double GetStepSizeNew() {return TIstepnew;} // necessary for time discrete input output elements 325 GetMinStepSize()326 virtual double GetMinStepSize() const {return TIminstep;} 327 328 GetStepRecommendation()329 virtual double GetStepRecommendation() const {return TIsteprecommendation;} SetStepRecommendation(double s)330 virtual void SetStepRecommendation(double s) {TIsteprecommendation = s;} 331 SetTime(double tt)332 virtual void SetTime(double tt) {TItime = tt; TIdrawtime = tt;} GetTime()333 virtual double GetTime() const {return TItime;} 334 GetStepEndTime()335 virtual double GetStepEndTime() const {return TItime+TIstep;} GetDrawTime()336 virtual double GetDrawTime() const {return TIdrawtime;} // for MBS interface GetActualDrawTime()337 virtual double GetActualDrawTime() const {return TIdrawtime;} // for WinCompDriverInterface SetDrawTime(double t)338 virtual void SetDrawTime(double t) {TIdrawtime = t;} GetDrawVector()339 virtual const Vector& GetDrawVector() const {return TIx0draw;} GetDrawVector()340 virtual Vector& GetDrawVector() {return TIx0draw;} GetDrawValue(int i)341 virtual const double& GetDrawValue(int i) const {return TIx0draw(i);} GetDrawValue(int i)342 virtual double& GetDrawValue(int i) {return TIx0draw(i);} 343 GetCompTime()344 virtual double GetCompTime() const {return TIcomptime;} 345 virtual double GetClockTime() const; GetStepChanges()346 virtual int GetStepChanges() const {return log.changestep;} GetJacCount()347 virtual int GetJacCount() const {return log.jaccount;} 348 349 //#write problem dependent solution data WriteSol()350 virtual void WriteSol() {} 351 352 //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 353 //main functions for data handling 354 355 // these functions are due to the data storage 356 // after each call of ComputationFeedBack::ResultsUpdated() 357 // the function StoreResultsIsOn() is called to check if the results should be stored 358 virtual int StoreResultsIsOn(); 359 360 // here the working object should remove all stored states 361 virtual int RemoveResults(); 362 363 // here the working object should save its actual state 364 virtual void StoreResults(DataSaver & storage, double & m_TimePoint); 365 366 // here the working object should replace its actual state with the one from the DataLoader 367 // the function is called from the driving module with the user request 368 virtual void LoadResults(DataLoader & loader, int m_TimePointNumber); 369 370 // this line identifies the current data storage structure 371 // when the data set or the way of storage is changed, 372 // new version identifier should be introduced 373 //virtual const char * GetDataFormatVersionIdentifier() {return "Version0.1";}; 374 375 //TIME-INT VARIABLES GetSolSet()376 virtual const SolverSettings& GetSolSet() const {return solset;} //$ YV 2012-11-29: MBSSolSet does not exist any longer GetSolSet()377 virtual SolverSettings& GetSolSet() {return solset;} SetSolSet(const SolverSettings & solsetI)378 virtual void SetSolSet(const SolverSettings& solsetI) {solset = solsetI;} 379 WriteSolTau()380 const double& WriteSolTau() const {return TIwritesolevery;} WriteSolTau()381 double& WriteSolTau() {return TIwritesolevery;} 382 MaxDiscontinuousIt()383 int MaxDiscontinuousIt() const {return solset.maxdiscontinuousit;}; MaxDiscontinuousIt()384 int& MaxDiscontinuousIt() {return solset.maxdiscontinuousit;}; 385 DiscontinuousAccuracy()386 double DiscontinuousAccuracy() const {return solset.discontinuousaccuracy;}; DiscontinuousAccuracy()387 double& DiscontinuousAccuracy() {return solset.discontinuousaccuracy;}; 388 389 //virtual int StopCalculation() const {return 1; }; StopCalculation()390 virtual int StopCalculation() const {return bStopWhenPossible;}; PauseCalculation()391 virtual int PauseCalculation() const {return bPause;}; 392 AbsAccuracy()393 double AbsAccuracy() const {return solset.absaccuracy;}; AbsAccuracy()394 double& AbsAccuracy() {return solset.absaccuracy;}; 395 FullAdaptive()396 int FullAdaptive() const {return solset.fulladaptive;}; FullAdaptive()397 int& FullAdaptive() {return solset.fulladaptive;}; 398 KVersion()399 int KVersion() const {return TIkv;}; KVersion()400 int& KVersion() {return TIkv;}; 401 SetStepDiscontinuous()402 void SetStepDiscontinuous() {TIdiscontstep++;} IsStepDiscontinuous()403 int IsStepDiscontinuous() const {return (TIdiscontstep!=0);} 404 IsJacobianComputation()405 virtual int IsJacobianComputation() const {return computejacflag;} SetJacobianComputationFlag(int i)406 virtual void SetJacobianComputationFlag(int i) {computejacflag = i;} 407 NLS_UseSparseSolver()408 virtual int NLS_UseSparseSolver() const {return solset.nls_usesparsesolver;}; NLS_UseSparseSolver()409 virtual int& NLS_UseSparseSolver() {return solset.nls_usesparsesolver;}; 410 NLS_SolveUndeterminedSystem()411 virtual int NLS_SolveUndeterminedSystem() const {return solset.nls_solve_undetermined_system;}; NLS_EstimatedConditionNumber()412 virtual double NLS_EstimatedConditionNumber() const {return solset.nls_estimated_condition_number;} 413 414 //NEWTON-METHOD: NLS_ModifiedNewton()415 int NLS_ModifiedNewton() const {return solset.nls_modifiednewton;}; NLS_ModifiedNewton()416 int& NLS_ModifiedNewton() {return solset.nls_modifiednewton;}; 417 NLS_AbsoluteAccuracy()418 double NLS_AbsoluteAccuracy() const {return solset.nls_absoluteaccuracy;}; NLS_AbsoluteAccuracy()419 double& NLS_AbsoluteAccuracy() {return solset.nls_absoluteaccuracy;}; 420 NLS_RelativeAccuracy()421 double NLS_RelativeAccuracy() const {return solset.nls_relativeaccuracy;}; NLS_RelativeAccuracy()422 double& NLS_RelativeAccuracy() {return solset.nls_relativeaccuracy;}; 423 NLS_NumDiffepsi()424 double NLS_NumDiffepsi() const {return solset.nls_numdiffepsi;}; NLS_NumDiffepsi()425 double& NLS_NumDiffepsi() {return solset.nls_numdiffepsi;}; 426 NLS_MaxNewtonSteps()427 int NLS_MaxNewtonSteps() const {return solset.nls_maxmodnewtonsteps + solset.nls_maxrestartnewtonsteps + solset.nls_maxfullnewtonsteps;}; 428 NLS_MaxModNewtonSteps()429 int NLS_MaxModNewtonSteps() const {return solset.nls_maxmodnewtonsteps;}; NLS_MaxModNewtonSteps()430 int& NLS_MaxModNewtonSteps() {return solset.nls_maxmodnewtonsteps;}; 431 NLS_MaxRestartNewtonSteps()432 int NLS_MaxRestartNewtonSteps() const {return solset.nls_maxrestartnewtonsteps;}; NLS_MaxRestartNewtonSteps()433 int& NLS_MaxRestartNewtonSteps() {return solset.nls_maxrestartnewtonsteps;}; 434 NLS_MaxFullNewtonSteps()435 int NLS_MaxFullNewtonSteps() const {return solset.nls_maxfullnewtonsteps;}; NLS_MaxFullNewtonSteps()436 int& NLS_MaxFullNewtonSteps() {return solset.nls_maxfullnewtonsteps;}; 437 NLS_TrustRegion()438 int NLS_TrustRegion() const {return solset.nls_trustregion;}; NLS_TrustRegion()439 int& NLS_TrustRegion() {return solset.nls_trustregion;}; 440 NLS_TrustRegionDiv()441 double NLS_TrustRegionDiv() const {return solset.nls_trustregiondiv;}; NLS_TrustRegionDiv()442 double& NLS_TrustRegionDiv() {return solset.nls_trustregiondiv;}; 443 NLS_SymmetricJacobian()444 int NLS_SymmetricJacobian() const {return solset.nls_symmetricjacobian;}; NLS_SymmetricJacobian()445 int& NLS_SymmetricJacobian() {return solset.nls_symmetricjacobian;}; 446 NumSolver()447 const NumNLSolver& NumSolver() const {return numsol;} NumSolver()448 NumNLSolver& NumSolver() {return numsol;} 449 //void RestartNumSolver() {numsol = NumNLSolver(this);} 450 451 //void SetJacCount() {jaccount++;} SetJacCount(int ii)452 void SetJacCount(int ii) {log.jaccount=ii;} 453 454 //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 455 //#graphics: DrawSystem()456 virtual void DrawSystem() {}; GetRC()457 virtual RenderContext * GetRC() { return pCurrentRC; } GetCWC()458 virtual ControlWindowContext * GetCWC() { return p2DDrawWindow; } 459 virtual void MyDrawLineH(const Vector3D& p1, const Vector3D&p2, const Vector3D& vy2, 460 double t, double h, int drawouterface=1) const; 461 462 void DrawLegend(double yoff=0) const; 463 464 void DrawColorQuads(const TArray<Vector3D>& p, const TArray<double>& v, int n1, int n2, 465 int colormode, int drawlines=0, int vres=1); //draw n1*n2 quads with FEcolors, interpolate values v with higher resultion vres 466 467 void DrawQuad(const Vector3D& p1,const Vector3D& p2,const Vector3D& p3,const Vector3D& p4) const; 468 469 void DrawTrig(const Vector3D& p1,const Vector3D& p2,const Vector3D& p3) const; 470 471 void DrawHex(const Vector3D& p1,const Vector3D& p2,const Vector3D& p3,const Vector3D& p4, 472 const Vector3D& p5,const Vector3D& p6,const Vector3D& p7,const Vector3D& p8, int drawouterfaces=1) const; 473 474 void DrawCube(const Vector3D& p0, const Vector3D& v1, const Vector3D& v2, const Vector3D& v3) const; //draws cube with reference point and 3 axes 475 476 477 virtual void MyDrawRectangle(const Vector3D& p1, const Vector3D& p2, const Vector3D& p3, const Vector3D& p4, double thickness, const Vector3D* colline, const Vector3D* colfill=0) const; 478 479 virtual void MyDrawLine(const Vector3D& p1, const Vector3D&p2, double thickness) const; 480 virtual void MyDrawLine(const Vector3D& p1, const Vector3D&p2, double thickness, const Vector3D& col) const; MyDrawLine(const Vector3D & p1,const Vector3D & p2,double t,double h)481 virtual void MyDrawLine(const Vector3D& p1, const Vector3D&p2, double t, double h) const 482 { 483 MyDrawLineH(p1,p2,p2-p1,t,h); 484 } 485 virtual void MyDrawCircleXY(const Vector3D p, double r, const Vector3D& col, int res=12, double thickness=1) const; 486 virtual void MyDrawArrow(const Vector3D& p1, const Vector3D&p2, const Vector3D& col, 487 double linethickness = -1, double headsize = -1, int resolution=8) const; 488 void DrawArrow(const Vector3D& p1, const Vector3D&p2, double linethickness = -1, double headsize = -1, int resolution = 8) const; 489 void DrawColorArrow(double v, const Vector3D& p1, const Vector3D&p2, double linethickness = -1, double headsize = -1, int resolution = 8); 490 491 492 virtual void DrawColorZyl(double v, const Vector3D& pz1,const Vector3D& pz2, double rz, int tile); //draw a zylinder in 3D, with color specified by value v, endpoints p1 and p2, Radius r, number of surfaces (discretized): tile 493 virtual void DrawZyl(const Vector3D& pz1, const Vector3D& pz2, double rz, int tile=8) const; 494 virtual void DrawCone(const Vector3D& pz1, const Vector3D& pz2, double rz, int tile=8, int drawconelines = 0) const; 495 virtual void DrawZyl(const Vector3D& pz1, const Vector3D& pz2, const Vector3D& pz1dir,const Vector3D& pz2dir, double rz, int leftend = 1, int rightend = 1, int tile=8) const; 496 497 //tile ... number of tiles per circumference, fill:1 ...draw whole sphere, fill=0.5 ... draw half sphere 498 virtual void DrawSphere(const Vector3D& p1, double r, int tile=8, double fill=1) const; 499 virtual void DrawColorSphere(double value, const Vector3D& p, double r, int tile=8, double fill=1); 500 501 502 virtual void DrawPolygon(const TArray<Vector3D>& p, int drawlines=0, double linewidth=1) const; 503 virtual void DrawPolygonOutline(const TArray<Vector3D>& p, double linewidth=1) const; 504 505 virtual void SetColor(const Vector3D& col); //set new color (set actcolor) SetLineColor(const Vector3D & col)506 virtual void SetLineColor(const Vector3D& col) {colline = col;}; //set new color (set actcolor) GetLineColor()507 virtual const Vector3D& GetLineColor() {return colline;}; //set new color (set actcolor) SetTransparency(int transp_on)508 virtual void SetTransparency(int transp_on) {transparency_on = transp_on;}; GetTransparency()509 virtual int GetTransparency() const {return transparency_on;}; 510 virtual void ChooseColor(float R, float G, float B) const; //choose color for painting SetDrawlines(int i)511 virtual void SetDrawlines(int i) {drawlines = i;} SetDrawlinesH(int i)512 virtual void SetDrawlinesH(int i) {drawlinesh = i;} 513 //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 514 //!AD: 2012-12-13: functions for 2D Draw window DrawSystem2D()515 virtual void DrawSystem2D() {}; 516 //virtual void AddDrawComponent(ControlWindowContext::DrawComponent dcomp); 517 virtual void AddDrawComponent_Line(int mbs_elnr, TMBSElementType mbs_type, int sub_elnr, TDrawElementType sub_type, Vector2D p1, Vector2D p2, Vector3D col); 518 virtual void AddDrawComponent_Rect(int mbs_elnr, TMBSElementType mbs_type, int sub_elnr, TDrawElementType sub_type, Vector2D center, Vector2D size, Vector3D col_border, Vector3D col_background); 519 virtual void AddDrawComponent_Ellipse(int mbs_elnr, TMBSElementType mbs_type, int sub_elnr, TDrawElementType sub_type, Vector2D center, Vector2D size, Vector3D col_border, Vector3D col_background); 520 virtual void AddDrawComponent_Text(int mbs_elnr, TMBSElementType mbs_type, int sub_elnr, TDrawElementType sub_type, Vector2D center, Vector2D size, Vector3D col_text, mystr& text, TTextAllign positioning); 521 // functions that allow the 2D Draw Window to change properties of MBS Elements 522 //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 523 GetDrawResolution()524 virtual int GetDrawResolution() const {return GetIOption(139);} //depreciated --> use options! 525 526 //to be overwritten in MBS: 527 //get values in MBS_EDC: MBS_EDC_TreeGetDouble(double & data,const char * name)528 virtual void MBS_EDC_TreeGetDouble(double& data, const char* name) const {}; MBS_EDC_TreeGetInt(int & data,const char * name)529 virtual void MBS_EDC_TreeGetInt(int& data, const char* name) const {}; MBS_EDC_TreeGetString(const char * name)530 virtual const char* MBS_EDC_TreeGetString(const char* name) const {return 0;}; 531 //set values in MBS_EDC: MBS_EDC_TreeSetDouble(double data,const char * name)532 virtual void MBS_EDC_TreeSetDouble(double data, const char* name) {}; MBS_EDC_TreeSetInt(int data,const char * name)533 virtual void MBS_EDC_TreeSetInt(int data, const char* name) {}; MBS_EDC_TreeSetString(char * data,const char * name)534 virtual void MBS_EDC_TreeSetString(char* data, const char* name) {}; 535 536 GetTImincol()537 virtual const double& GetTImincol() const {return TImincol;} GetFEmincol()538 virtual double& GetFEmincol() {return TImincol;} GetTImaxcol()539 virtual const double& GetTImaxcol() const {return TImaxcol;} GetFEmaxcol()540 virtual double& GetFEmaxcol() {return TImaxcol;} IsFlagComputeMinMaxFECol()541 virtual int IsFlagComputeMinMaxFECol() const {return flag_compute_minmax_FEcol;} SetFlagComputeMinMaxFECol(int flag)542 virtual void SetFlagComputeMinMaxFECol(int flag) {flag_compute_minmax_FEcol = flag;} 543 544 //+++++++++++++++++++++++++++++++++++++++++++++++++ 545 //element data access: GetNElements()546 virtual int GetNElements() const {return 0;} //overwritten in mbs.h GetElementData(int i,int type,int value,ElementDataContainer & edc)547 virtual void GetElementData(int i, int type, int value, ElementDataContainer& edc) {}; SetElementData(int i,int type,int value,ElementDataContainer & edc)548 virtual int SetElementData(int i, int type, int value, ElementDataContainer& edc) {return 0;}; 549 550 //$AD 2013-07-08: Manipulate Content of arrays in IOElements from 2D Draw Window 551 virtual void InsertIOElemConNode(int elemnr, int list_idx, int input_nr, Vector2D pos) = 0; 552 virtual void DeleteIOElemConNode(int elemnr, int list_idx) = 0; 553 554 //function call: 555 //virtual int CallCompFunction(int action, int option, int value, ElementDataContainer* edc) {return 0;}; 556 virtual int CallCompFunction(int action, int option = 0, int value = 0, ElementDataContainer* edc = NULL) {return 0;}; 557 558 //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 559 // access functions for Sensor / PlotTool (Plottool has no classes MBS, Sensor, etc included) 560 virtual TArray<double>* GetSensorValuesArrayPtr(int i) =0; // get the i-th sensors (as registered in mbs) own stored values from internal array 561 virtual TArray<double>* GetSensorTimesArrayPtr(int i) =0; // get the i-th sensors (as registered in mbs) own stored times from internal array GetSensorName(int i)562 virtual mystr GetSensorName(int i) {return 0;}; // get the i-th sensors (as registered in mbs) name GetSensorNumber(mystr & name)563 virtual int GetSensorNumber(mystr& name) {return 0;}; // get the number of the (first) sensor that matches the sensorname 564 565 //+++++++++++++++++++++++++++++++++++++++++++++++++ 566 virtual Vector3D FEColor(double val) const; 567 virtual void SetFEColor(double val) const; 568 virtual void DrawIsoQuad(Vector3D* p, const Vector3D& n, double* v); //draw quad with iso-lines 569 virtual void DrawIsoQuadTex(Vector3D* p, const Vector3D& n, double* v, int res, char* teximage); 570 571 // update minimum and maximum values by val 572 virtual void UpdateFEMinMaxCol(double val); 573 574 virtual void SetTexStoreResolution(int n); GetTexStoreResolution()575 virtual int GetTexStoreResolution() const {return texstorenn;} 576 virtual int GetTexStoreMaxResolution(); 577 578 //drawing: 579 // int TIdrawmode; //*YV: commented out as it is not needed any longer 580 Vector3D actcolor, colline; ColLine()581 virtual Vector3D & ColLine() { return colline; } // for interface access 582 int drawlines; 583 int drawlinesh; 584 int transparency_on; //activate / deactivate transpareny mode, factor is an option 585 586 //texture: 587 int TIdrawtexres; 588 int texstorenn; //size of texture, minimum requirement: 64x64 for opengl 589 int reduceimage; //reduce size of image, such that only a part of it is used 590 char* texImage; //pointer to texture image 591 unsigned int texName; 592 593 //public Time-Integration variables: 594 double ada_err_sum; 595 double acterror, acterror_damp; 596 double comperr; 597 double TItimeperstep; 598 int TIissubstep; //in fulladaptive, currently computing the 599 double oldenergy; 600 double TIcondnum; //rough estimate condition number of jacobian 601 602 //drawing/storing/writing results: 603 double laststoredata; 604 double lastloadresults; 605 int drawnow; 606 607 int system_matrices_written; //flag, system matrices are written at first time step 608 609 //graphics: 610 RenderContext * pCurrentRC; 611 ControlWindowContext* p2DDrawWindow; //!AD: new 2012-12-13 612 613 double tidraw_offx; 614 double tidraw_offy; 615 double tidraw_offz; 616 617 618 //PG: 619 //// output in log file 620 //virtual ofstream& LogOut() { return logout; } 621 ofstream logout; // log-file 622 mystr logout_name; // log-file name 623 624 protected: 625 TimeIntLog log; 626 627 protected: 628 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 629 //solver settings, can be accessed in derived class: 630 SolverSettings solset; 631 632 private: 633 RenderScene(RenderContext * pRC)634 virtual void RenderScene(RenderContext * pRC) 635 { 636 pCurrentRC = pRC; 637 if (!solset.withgraphics){return;}; 638 639 DrawSystem(); 640 }; 641 RenderControlWindow(ControlWindowContext * pCWC)642 virtual void RenderControlWindow(ControlWindowContext* pCWC) 643 { 644 p2DDrawWindow = pCWC ; //!AD: new 2012-12-13 645 if (!solset.withgraphics){return;}; 646 647 DrawSystem2D(); 648 } 649 650 651 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 652 653 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 654 //new structure vor TimeInt variables during computation 655 TimeIntVars TIV; 656 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 657 658 659 protected: 660 //erase: 661 //Vector TIx0last, TIk0last; 662 //Vector TINLsubstep_k0save, TINLsubstep_x0save; 663 //TIx0old, TIx0save, TIk0old, TIk0save; 664 //double TItimesave; //actual simulation time 665 666 //general time ingetration: 667 //actual state, old state (NonlinSubStep), saved state (Save/RestoreState), state for drawing 668 Vector TIx0; //actual states for computation: [q v x z], lengths [sos sos es is] 669 private: 670 Vector TIk0; //actual states for k-version of time integration 671 Vector TIdata; //data variables for actual step 672 double TItime; //actual time for computation 673 674 Vector TIx0draw; //actual states for drawing; 675 Vector TIdrawdata; //data variables for drawing; 676 double TIdrawtime; //actual time for drawing result; 677 678 Vector TIx0start; //starting vector for time integration and static solver 679 Vector TIk0start; //k-version of starting vector for time integration and static solver 680 681 TIstepsave TIlaststep_state; //state at end of last successful time step 682 TIstepsave TIlastnonlinit_state;//state after last nonlinear iteration 683 TIstepsave TItemp_state; //temporary state (for postprocessing) // AP: AND uses this state for the plastic strains after transport, before nonlinear iteration 684 685 double TIstep; //actual timestepsize 686 double TIstepnew; //suggested timestepsize for next timestep 687 double TImaxstep; //max timestepsize 688 double TIminstep; //min timestepsize 689 690 691 double TIsteprecommendation; //step suggested by mbs system 692 693 protected: 694 double TIcomptime; //elapsed computational time 695 double TIactlocerr, TInonls; 696 int TIit; 697 int TIrejectedsteps; 698 int TIdiscontstep; 699 700 private: 701 int TIfinished; //can be set 1 if integration shall be stopped; set -1 if integration stopped with error 702 703 double TImincol; 704 double TImaxcol; 705 int flag_compute_minmax_FEcol; 706 707 708 double loadfact; //for static computation, do not apply all load at once ..., for dynamic computation == 1 !!!! 709 int computejacflag; //is 1 during computation of jacobian 710 711 //Tableau data, variable stepsize 712 TArray<int> stage_tableaus; 713 TArray<IRK_Tableau*> tableaus; 714 int act_tableau; //actual used table 715 int act_stage; 716 717 //variables to control order selection 718 int steps_since_orderchange; 719 int min_newtonits; 720 int newjacobians_since_orderchange; 721 int last_jaccount; 722 int reduce_order; //reduce as soon as possible 723 724 double TIlastwritesol; 725 double TIwritesolevery; 726 int TInumber_of_solution_data_steps; //$ PG 2012-4-13: store number of stored solution data steps 727 728 //old: 729 int TInoincs; 730 double TIbasestep; 731 int impl_equ_flag; //for implicit differential equations 732 int nls_output; 733 734 protected: 735 NumNLSolver numsol; 736 737 private: 738 int TIkv; //K-version integration 739 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 740 //local vectors/matrices for Jacobian/LocalJacobian 741 SparseMatrix mlocgs; 742 SparseMatrix mlocf2s; 743 SparseMatrix mlocms; 744 745 Matrix mlocf; 746 Matrix mlocg; 747 Matrix mlocf2; 748 Matrix mlocm; 749 Matrix mlocm2; 750 Matrix diag; 751 Matrix TIm; //for K-version NLF 752 SparseMatrix TIm_sparse; 753 754 Vector locjacf0,locjacf20,locjacg0; //for non-symmetric jacobian 755 Vector locjacf1,locjacf21,locjacg1; 756 Vector locjacf2,locjacf22,locjacg2, xx; 757 SparseVector slocjacf20, slocjacf21, slocjacf22; 758 SparseVector slocjacg0, slocjacg1, slocjacg2; 759 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 760 761 int bandwidthm; 762 int reducedbandsize; //reduce bandsize-part of the matrix for certain flexible ACRS models connected with rigid bodies 763 int TIm_initialized; //initialize only at beginning 764 765 IVector temprowind, tempclearind; 766 //local vectors/matrices for GeneralIstep 767 768 //local vectors/matrices for NLF 769 Vector fg; 770 Vector gv[global_maxstages]; 771 Vector xv[global_maxstages]; 772 Vector iv[global_maxstages]; 773 Vector gu[global_maxstages]; 774 Vector gue[global_maxstages]; 775 Vector xi[global_maxstages]; 776 Vector kue[global_maxstages]; //first order 777 Vector ku[global_maxstages]; //second order u 778 Vector kv[global_maxstages]; //second order u 779 Vector xh; 780 Vector xh2; 781 Vector evalfs[global_maxstages]; 782 Vector evalf; 783 Vector evalfe; 784 Vector evalg; 785 Vector TIu0; 786 Vector TIv0; 787 Vector TIu0e; 788 Vector TIi0; 789 Vector TIku0; 790 Vector TIku0e; 791 Vector TIkv0; 792 Vector minvf[global_maxstages]; 793 //Matrix m[global_maxstages]; 794 }; 795 796 #endif 797