Ice Sheet System Model  4.18
Code documentation
Functions
SystemMatricesx.cpp File Reference
#include "./SystemMatricesx.h"
#include "../../shared/shared.h"
#include "../../toolkits/toolkits.h"
#include "../AllocateSystemMatricesx/AllocateSystemMatricesx.h"

Go to the source code of this file.

Functions

void SystemMatricesx (Matrix< IssmDouble > **pKff, Matrix< IssmDouble > **pKfs, Vector< IssmDouble > **ppf, Vector< IssmDouble > **pdf, IssmDouble *pkmax, FemModel *femmodel, bool isAllocated)
 

Function Documentation

◆ SystemMatricesx()

void SystemMatricesx ( Matrix< IssmDouble > **  pKff,
Matrix< IssmDouble > **  pKfs,
Vector< IssmDouble > **  ppf,
Vector< IssmDouble > **  pdf,
IssmDouble pkmax,
FemModel femmodel,
bool  isAllocated 
)

Definition at line 10 of file SystemMatricesx.cpp.

10  {
11 
12  /*intermediary: */
13  int i,M,N;
14  int analysisenum;
15  Element *element = NULL;
16  Load *load = NULL;
17 
18  /*output: */
19  Matrix<IssmDouble> *Kff = NULL;
20  Matrix<IssmDouble> *Kfs = NULL;
21  Vector<IssmDouble> *pf = NULL;
22  Vector<IssmDouble> *df = NULL;
23  IssmDouble kmax = 0;
24 
25  /*Display message*/
26  if(VerboseModule()) _printf0_(" Allocating matrices");
27 
28  /*retrieve parameters: */
30  Analysis* analysis = EnumToAnalysis(analysisenum);
31 
32  /*Check if there are penalties*/
33  bool ispenalty = femmodel->loads->IsPenalty();
34 
35  /*First, we might need to do a dry run to get kmax if penalties are employed*/
36  if(ispenalty){
37 
38  /*Allocate Kff_temp*/
39  Matrix<IssmDouble> *Kff_temp = NULL;
40  AllocateSystemMatricesx(&Kff_temp,NULL,NULL,NULL,femmodel);
41 
42  /*Get complete stiffness matrix without penalties*/
43  for (i=0;i<femmodel->elements->Size();i++){
44  element=xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i));
45  if(!element->AnyFSet() && analysisenum!=StressbalanceAnalysisEnum) continue;
46  ElementMatrix* Ke = analysis->CreateKMatrix(element);
47  ElementVector* pe = analysis->CreatePVector(element);
48  element->ReduceMatrices(Ke,pe);
49  if(Ke) Ke->AddToGlobal(Kff_temp,NULL);
50  delete Ke;
51  delete pe;
52  }
53 
54  for (i=0;i<femmodel->loads->Size();i++){
55  load=xDynamicCast<Load*>(femmodel->loads->GetObjectByOffset(i));
56  load->CreateKMatrix(Kff_temp,NULL);
57  }
58  Kff_temp->Assemble();
59 
60  /*Now, figure out maximum value of stiffness matrix, so that we can penalize it correctly: */
61  kmax=Kff_temp->Norm(NORM_INF);
62  delete Kff_temp;
63  }
64 
65  /*Allocate stiffness matrices and load vector*/
66  if(isAllocated) {
67  Kff = *pKff;
68  Kfs = *pKfs;
69  pf = *ppf;
70  df = *pdf;
71  }
72  else {
73  AllocateSystemMatricesx(&Kff,&Kfs,&df,&pf,femmodel);
74  }
75 
76  /*Display size*/
77  if(VerboseModule()){
78  Kff->GetSize(&M,&N);
79  _printf0_(" (Kff stiffness matrix size: "<<M<<" x "<<N<<")\n");
80  }
81 
82  if(VerboseModule()) _printf0_(" Assembling matrices\n");
83 
84  /*Fill stiffness matrix and load vector from elements*/
85  for (i=0;i<femmodel->elements->Size();i++){
86  element=xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i));
87  if(!element->AnyFSet() && analysisenum!=StressbalanceAnalysisEnum) continue;
88  ElementMatrix* Ke = analysis->CreateKMatrix(element);
89  ElementVector* pe = analysis->CreatePVector(element);
90  element->ReduceMatrices(Ke,pe);
91  if(Ke) Ke->AddToGlobal(Kff,Kfs);
92  if(pe){
93  pe->AddToGlobal(pf);
94  }
95  delete Ke;
96  delete pe;
97  }
98 
99  /*Fill stiffness matrix and load vector from loads*/
100  for(i=0;i<femmodel->loads->Size();i++){
101  load=xDynamicCast<Load*>(femmodel->loads->GetObjectByOffset(i));
102  load->CreateKMatrix(Kff,Kfs);
103  load->CreatePVector(pf);
104  }
105 
106  /*Now deal with penalties (only in loads)*/
107  if(ispenalty){
108  for (i=0;i<femmodel->loads->Size();i++){
109  load=xDynamicCast<Load*>(femmodel->loads->GetObjectByOffset(i));
110  load->PenaltyCreateKMatrix(Kff,Kfs,kmax);
111  load->PenaltyCreatePVector(pf,kmax);
112  }
113  }
114 
115  /*Create dof vector for stiffness matrix preconditioning*/
116  if(pdf){
117  for(i=0;i<femmodel->elements->Size();i++){
118  element=xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i));
119  ElementVector* de=analysis->CreateDVector(element);
120  if(de) de->InsertIntoGlobal(df);
121  delete de;
122  }
123  }
124 
125  /*Assemble matrices and vector*/
126  Kff->Assemble();
127  Kfs->Assemble();
128  pf->Assemble();
129  df->Assemble();
130  //Kff->AllocationInfo();
131  //Kfs->AllocationInfo();
132 
133  /*cleanu up and assign output pointers: */
134  delete analysis;
135  if(pKff) *pKff=Kff;
136  else delete Kff;
137  if(pKfs) *pKfs=Kfs;
138  else delete Kfs;
139  if(ppf) *ppf=pf;
140  else delete pf;
141  if(pdf) *pdf=df;
142  else delete df;
143  if(pkmax) *pkmax=kmax;
144 }
DataSet::Size
int Size()
Definition: DataSet.cpp:399
Matrix< IssmDouble >
IssmDouble
double IssmDouble
Definition: types.h:37
_printf0_
#define _printf0_(StreamArgs)
Definition: Print.h:29
StressbalanceAnalysisEnum
@ StressbalanceAnalysisEnum
Definition: EnumDefinitions.h:1285
Load::CreateKMatrix
virtual void CreateKMatrix(Matrix< IssmDouble > *Kff, Matrix< IssmDouble > *Kfs)=0
FemModel::parameters
Parameters * parameters
Definition: FemModel.h:46
Matrix::Assemble
void Assemble(void)
Definition: Matrix.h:185
NORM_INF
@ NORM_INF
Definition: toolkitsenums.h:15
Load
Definition: Load.h:22
AllocateSystemMatricesx
void AllocateSystemMatricesx(Matrix< IssmDouble > **pKff, Matrix< IssmDouble > **pKfs, Vector< IssmDouble > **pdf, Vector< IssmDouble > **ppf, FemModel *femmodel)
Definition: AllocateSystemMatricesx.cpp:9
VerboseModule
bool VerboseModule(void)
Definition: Verbosity.cpp:23
Analysis::CreatePVector
virtual ElementVector * CreatePVector(Element *element)=0
Element
Definition: Element.h:41
ElementVector::InsertIntoGlobal
void InsertIntoGlobal(Vector< IssmDouble > *pf)
Definition: ElementVector.cpp:238
Load::PenaltyCreatePVector
virtual void PenaltyCreatePVector(Vector< IssmDouble > *pf, IssmDouble kmax)=0
Element::ReduceMatrices
virtual void ReduceMatrices(ElementMatrix *Ke, ElementVector *pe)=0
Loads::IsPenalty
bool IsPenalty()
Definition: Loads.cpp:124
FemModel::loads
Loads * loads
Definition: FemModel.h:54
FemModel::elements
Elements * elements
Definition: FemModel.h:44
Load::CreatePVector
virtual void CreatePVector(Vector< IssmDouble > *pf)=0
Vector::Assemble
void Assemble(void)
Definition: Vector.h:142
EnumToAnalysis
Analysis * EnumToAnalysis(int analysis_enum)
Definition: EnumToAnalysis.cpp:13
Matrix::Norm
IssmDouble Norm(NormMode norm_type)
Definition: Matrix.h:197
Element::AnyFSet
bool AnyFSet(void)
Definition: Element.cpp:51
DataSet::GetObjectByOffset
Object * GetObjectByOffset(int offset)
Definition: DataSet.cpp:334
Analysis::CreateDVector
virtual ElementVector * CreateDVector(Element *element)=0
Parameters::FindParam
void FindParam(bool *pinteger, int enum_type)
Definition: Parameters.cpp:262
ElementVector::AddToGlobal
void AddToGlobal(Vector< IssmDouble > *pf)
Definition: ElementVector.cpp:155
Load::PenaltyCreateKMatrix
virtual void PenaltyCreateKMatrix(Matrix< IssmDouble > *Kff, Matrix< IssmDouble > *Kfs, IssmDouble kmax)=0
AnalysisTypeEnum
@ AnalysisTypeEnum
Definition: EnumDefinitions.h:36
ElementMatrix::AddToGlobal
void AddToGlobal(Matrix< IssmDouble > *Kff, Matrix< IssmDouble > *Kfs)
Definition: ElementMatrix.cpp:271
ElementVector
Definition: ElementVector.h:20
ElementMatrix
Definition: ElementMatrix.h:19
Vector< IssmDouble >
Analysis
Definition: Analysis.h:30
Matrix::GetSize
void GetSize(int *pM, int *pN)
Definition: Matrix.h:213
Analysis::CreateKMatrix
virtual ElementMatrix * CreateKMatrix(Element *element)=0
femmodel
FemModel * femmodel
Definition: esmfbinders.cpp:16