 |
Ice Sheet System Model
4.18
Code documentation
|
#include <StressbalanceAnalysis.h>
|
void | CreateConstraints (Constraints *constraints, IoModel *iomodel) |
|
void | CreateLoads (Loads *loads, IoModel *iomodel) |
|
void | CreateNodes (Nodes *nodes, IoModel *iomodel, bool isamr=false) |
|
int | DofsPerNode (int **doflist, int domaintype, int approximation) |
|
void | UpdateElements (Elements *elements, Inputs2 *inputs2, IoModel *iomodel, int analysis_counter, int analysis_type) |
|
void | UpdateParameters (Parameters *parameters, IoModel *iomodel, int solution_enum, int analysis_enum) |
|
void | Core (FemModel *femmodel) |
|
ElementVector * | CreateDVector (Element *element) |
|
ElementMatrix * | CreateJacobianMatrix (Element *element) |
|
ElementMatrix * | CreateKMatrix (Element *element) |
|
ElementVector * | CreatePVector (Element *element) |
|
void | GetSolutionFromInputs (Vector< IssmDouble > *solution, Element *element) |
|
void | GetSolutionFromInputsHoriz (Vector< IssmDouble > *solution, Element *element) |
|
void | GradientJ (Vector< IssmDouble > *gradient, Element *element, int control_type, int control_index) |
|
void | InputUpdateFromSolution (IssmDouble *solution, Element *element) |
|
void | UpdateConstraints (FemModel *femmodel) |
|
ElementMatrix * | CreateJacobianMatrixSSA (Element *element) |
|
ElementMatrix * | CreateKMatrixSSA (Element *element) |
|
ElementMatrix * | CreateKMatrixSSAFriction (Element *element) |
|
ElementMatrix * | CreateKMatrixSSALateralFriction (Element *element) |
|
ElementMatrix * | CreateKMatrixSSAViscous (Element *element) |
|
ElementVector * | CreatePVectorSSA (Element *element) |
|
ElementVector * | CreatePVectorSSAFront (Element *element) |
|
ElementVector * | CreatePVectorSSADrivingStress (Element *element) |
|
void | GetBSSA (IssmDouble *B, Element *element, int dim, IssmDouble *xyz_list, Gauss *gauss) |
|
void | GetBSSAFriction (IssmDouble *B, Element *element, int dim, IssmDouble *xyz_list, Gauss *gauss) |
|
void | GetBSSAprime (IssmDouble *B, Element *element, int dim, IssmDouble *xyz_list, Gauss *gauss) |
|
void | InputUpdateFromSolutionSSA (IssmDouble *solution, Element *element) |
|
ElementMatrix * | CreateKMatrixL1L2 (Element *element) |
|
ElementMatrix * | CreateKMatrixL1L2Friction (Element *element) |
|
ElementMatrix * | CreateKMatrixL1L2Viscous (Element *element) |
|
ElementVector * | CreatePVectorL1L2 (Element *element) |
|
ElementVector * | CreatePVectorL1L2Front (Element *element) |
|
ElementVector * | CreatePVectorL1L2DrivingStress (Element *element) |
|
void | InputUpdateFromSolutionL1L2 (IssmDouble *solution, Element *element) |
|
ElementMatrix * | CreateJacobianMatrixHO (Element *element) |
|
ElementMatrix * | CreateKMatrixHO (Element *element) |
|
ElementMatrix * | CreateKMatrixHOFriction (Element *element) |
|
ElementMatrix * | CreateKMatrixHOViscous (Element *element) |
|
ElementVector * | CreatePVectorHO (Element *element) |
|
ElementVector * | CreatePVectorHOFront (Element *element) |
|
ElementVector * | CreatePVectorHODrivingStress (Element *element) |
|
void | GetBHO (IssmDouble *B, Element *element, int dim, IssmDouble *xyz_list, Gauss *gauss) |
|
void | GetBHOFriction (IssmDouble *B, Element *element, int dim, IssmDouble *xyz_list, Gauss *gauss) |
|
void | GetBHOprime (IssmDouble *B, Element *element, int dim, IssmDouble *xyz_list, Gauss *gauss) |
|
void | InputUpdateFromSolutionHO (IssmDouble *solution, Element *element) |
|
ElementVector * | CreateDVectorFS (Element *element) |
|
ElementMatrix * | CreateJacobianMatrixFS (Element *element) |
|
ElementMatrix * | CreateKMatrixFS (Element *element) |
|
ElementMatrix * | CreateKMatrixFSFriction (Element *element) |
|
ElementMatrix * | CreateKMatrixFSShelf (Element *element) |
|
ElementMatrix * | CreateKMatrixFSViscous (Element *element) |
|
ElementMatrix * | CreateKMatrixFSViscousLA (Element *element) |
|
ElementMatrix * | CreateKMatrixFSViscousXTH (Element *element) |
|
ElementMatrix * | CreatePressureMassMatrix (Element *element) |
|
ElementMatrix * | CreateSchurPrecondMatrix (Element *element) |
|
ElementVector * | CreatePVectorFS (Element *element) |
|
ElementVector * | CreatePVectorFSFriction (Element *element) |
|
ElementVector * | CreatePVectorFSFront (Element *element) |
|
ElementVector * | CreatePVectorFSShelf (Element *element) |
|
ElementVector * | CreatePVectorFSStress (Element *element) |
|
ElementVector * | CreatePVectorFSViscous (Element *element) |
|
ElementVector * | CreatePVectorFSViscousLA (Element *element) |
|
ElementVector * | CreatePVectorFSViscousXTH (Element *element) |
|
void | GetBFS (IssmDouble *B, Element *element, int dim, IssmDouble *xyz_list, Gauss *gauss) |
|
void | GetBFSFriction (IssmDouble *B, Element *element, int dim, IssmDouble *xyz_list, Gauss *gauss) |
|
void | GetBFSprime (IssmDouble *B, Element *element, int dim, IssmDouble *xyz_list, Gauss *gauss) |
|
void | GetBFSprimeUzawa (IssmDouble *B, Element *element, int dim, IssmDouble *xyz_list, Gauss *gauss) |
|
void | GetBFSprimevel (IssmDouble *B, Element *element, int dim, IssmDouble *xyz_list, Gauss *gauss) |
|
void | GetBFSUzawa (IssmDouble *B, Element *element, int dim, IssmDouble *xyz_list, Gauss *gauss) |
|
void | GetBFSvel (IssmDouble *B, Element *element, int dim, IssmDouble *xyz_list, Gauss *gauss) |
|
void | GetCFS (IssmDouble *B, Element *element, int dim, IssmDouble *xyz_list, Gauss *gauss) |
|
void | GetCFSprime (IssmDouble *B, Element *element, int dim, IssmDouble *xyz_list, Gauss *gauss) |
|
void | GetSolutionFromInputsFS (Vector< IssmDouble > *solution, Element *element) |
|
void | InitializeXTH (Elements *elements, Parameters *parameters) |
|
void | InputUpdateFromSolutionFS (IssmDouble *solution, Element *element) |
|
void | InputUpdateFromSolutionFSXTH_d (Elements *elements, Parameters *parameters) |
|
void | InputUpdateFromSolutionFSXTH_tau (Elements *elements, Parameters *parameters) |
|
ElementMatrix * | CreateKMatrixCouplingHOFS (Element *element) |
|
ElementMatrix * | CreateKMatrixCouplingSSAFS (Element *element) |
|
ElementMatrix * | CreateKMatrixCouplingSSAFSFriction (Element *element) |
|
ElementMatrix * | CreateKMatrixCouplingSSAFSViscous (Element *element) |
|
ElementMatrix * | CreateKMatrixCouplingSSAHO (Element *element) |
|
ElementMatrix * | CreateKMatrixCouplingSSAHOFriction (Element *element) |
|
ElementMatrix * | CreateKMatrixCouplingSSAHOViscous (Element *element) |
|
ElementMatrix * | CreateKMatrixHOFS (Element *element) |
|
ElementMatrix * | CreateKMatrixSSAFS (Element *element) |
|
ElementMatrix * | CreateKMatrixSSAHO (Element *element) |
|
ElementMatrix * | CreateKMatrixSSA3d (Element *element) |
|
ElementMatrix * | CreateKMatrixSSA3dFriction (Element *element) |
|
ElementMatrix * | CreateKMatrixSSA3dViscous (Element *element) |
|
ElementVector * | CreatePVectorSSAFS (Element *element) |
|
ElementVector * | CreatePVectorSSAHO (Element *element) |
|
ElementVector * | CreatePVectorCouplingSSAFS (Element *element) |
|
ElementVector * | CreatePVectorCouplingSSAFSFriction (Element *element) |
|
ElementVector * | CreatePVectorCouplingSSAFSViscous (Element *element) |
|
ElementVector * | CreatePVectorHOFS (Element *element) |
|
ElementVector * | CreatePVectorCouplingHOFS (Element *element) |
|
ElementVector * | CreatePVectorCouplingHOFSFriction (Element *element) |
|
ElementVector * | CreatePVectorCouplingHOFSViscous (Element *element) |
|
void | GetBprimeSSAFS (IssmDouble *Bprime, Element *element, IssmDouble *xyz_list, Gauss *gauss) |
|
void | GetBprimeSSAFSTria (IssmDouble *Bprime, Element *element, IssmDouble *xyz_list, Gauss *gauss) |
|
void | GetBSSAFS (IssmDouble *B, Element *element, IssmDouble *xyz_list, Gauss *gauss) |
|
void | GetBSSAFSTria (IssmDouble *B, Element *element, IssmDouble *xyz_list, Gauss *gauss) |
|
void | GetBSSAHO (IssmDouble *B, Element *element, IssmDouble *xyz_list, Gauss *gauss) |
|
void | GetLFSSSA (IssmDouble *L, Element *element, Gauss *gauss) |
|
void | GetLprimeFSSSA (IssmDouble *Lprime, Element *element, IssmDouble *xyz_list, Gauss *gauss) |
|
void | GetLprimeSSAFS (IssmDouble *Lprime, Element *element, IssmDouble *xyz_list, Gauss *gauss) |
|
void | GetLSSAFS (IssmDouble *L, Element *element, Gauss *gauss) |
|
void | InputUpdateFromSolutionHOFS (IssmDouble *solution, Element *element) |
|
void | InputUpdateFromSolutionSSAFS (IssmDouble *solution, Element *element) |
|
void | InputUpdateFromSolutionSSAHO (IssmDouble *solution, Element *element) |
|
virtual | ~Analysis () |
|
Definition at line 11 of file StressbalanceAnalysis.h.
◆ CreateConstraints()
void StressbalanceAnalysis::CreateConstraints |
( |
Constraints * |
constraints, |
|
|
IoModel * |
iomodel |
|
) |
| |
|
virtual |
Implements Analysis.
Definition at line 24 of file StressbalanceAnalysis.cpp.
28 int finiteelement,p_fe,v_fe;
32 bool isSIA,isSSA,isL1L2,isHO,isFS,iscoupling;
33 bool spcpresent =
false;
56 iomodel->
FindConstant(&FSreconditioning,
"md.stressbalance.FSreconditioning");
64 if(!isSSA && !isHO && !isFS && !isL1L2)
return;
67 if((isSIA?1.:0.) + (isSSA?1.:0.) + (isL1L2?1.:0.) + (isHO?1.:0.) + (isFS?1.:0.) >1.)
76 if(isSSA) iomodel->
FindConstant(&finiteelement,
"md.flowequation.fe_SSA");
77 else if(isL1L2) finiteelement =
P1Enum;
78 else if(isHO) iomodel->
FindConstant(&finiteelement,
"md.flowequation.fe_HO");
79 else if(isFS){ iomodel->
FindConstant(&finiteelement,
"md.flowequation.fe_FS");
81 switch(finiteelement){
96 _error_(
"model not supported yet");
102 iomodel->
FetchData(&vertices_type,NULL,NULL,
"md.flowequation.vertex_equation");
103 iomodel->
FetchData(&nodeonFS,NULL,NULL,
"md.flowequation.borderFS");
104 iomodel->
FetchData(&nodeonbase,NULL,NULL,
"md.mesh.vertexonbase");
105 iomodel->
FetchData(&groundedice_ls,NULL,NULL,
"md.mask.ocean_levelset");
107 iomodel->
FetchData(&spcvz,&Mz,&Nz,
"md.stressbalance.spcvz");
110 iomodel->
FetchData(&spcvz,&Mz,&Nz,
"md.stressbalance.spcvy");
119 iomodel->
DeleteData(spcvz,
"md.stressbalance.spcvz");
124 iomodel->
DeleteData(spcvz,
"md.stressbalance.spcvy");
129 iomodel->
DeleteData(vertices_type,
"md.flowequation.vertex_equation");
130 iomodel->
DeleteData(nodeonFS,
"md.flowequation.borderFS");
131 iomodel->
DeleteData(nodeonbase,
"md.mesh.vertexonbase");
132 iomodel->
DeleteData(groundedice_ls,
"md.mask.ocean_levelset");
135 int count = constraints->
Size();
136 iomodel->
FetchData(&vertices_type,NULL,NULL,
"md.flowequation.vertex_equation");
137 iomodel->
FetchData(&surface,NULL,NULL,
"md.geometry.surface");
138 iomodel->
FetchData(&z,NULL,NULL,
"md.mesh.z");
139 bool addpressurespc =
false;
142 addpressurespc =
true;
147 switch(finiteelement){
179 _error_(
"not implemented yet");
182 iomodel->
DeleteData(vertices_type,
"md.flowequation.vertex_equation");
183 iomodel->
DeleteData(surface,
"md.geometry.surface");
197 iomodel->
FetchData(&spcvx,&Mx,&Nx,
"md.stressbalance.spcvx");
198 iomodel->
FetchData(&spcvy,&My,&Ny,
"md.stressbalance.spcvy");
199 iomodel->
FetchData(&spcvz,&Mz,&Nz,
"md.stressbalance.spcvz");
200 iomodel->
FetchData(&nodeonSSA,NULL,NULL,
"md.flowequation.borderSSA");
205 iomodel->
FetchData(&vertices_type,NULL,NULL,
"md.flowequation.vertex_equation");
206 iomodel->
FetchData(&surface,NULL,NULL,
"md.geometry.surface");
207 iomodel->
FetchData(&z,NULL,NULL,
"md.mesh.z");
213 timesx=xNew<IssmDouble>(Nx);
215 timesx[j]=spcvx[(Mx-1)*Nx+j];
218 timesy=xNew<IssmDouble>(Ny);
220 timesy[j]=spcvy[(My-1)*Ny+j];
223 timesz=xNew<IssmDouble>(Nz);
225 timesz[j]=spcvz[(Mz-1)*Nz+j];
235 if (reCast<int,IssmDouble>(nodeonHO[i])){
240 if (!xIsNan<IssmDouble>(spcvx[i])){
244 if (!xIsNan<IssmDouble>(spcvy[i])){
250 else if (reCast<int,IssmDouble>(nodeonSSA[i])){
255 if (!xIsNan<IssmDouble>(spcvx[i])){
259 if (!xIsNan<IssmDouble>(spcvy[i])){
265 else _error_(
"if vertices_type is SSAHO, you shoud have nodeonHO or nodeonSSA");
270 if (reCast<int,IssmDouble>(nodeonHO[i])){
277 if (!xIsNan<IssmDouble>(spcvx[i])){
281 if (!xIsNan<IssmDouble>(spcvy[i])){
287 else if (reCast<int,IssmDouble>(nodeonFS[i])){
292 if (!xIsNan<IssmDouble>(spcvx[i])){
296 if (!xIsNan<IssmDouble>(spcvy[i])){
300 if (!xIsNan<IssmDouble>(spcvz[i])){
305 else _error_(
"if vertices_type is HOFS, you shoud have nodeonHO or nodeonFS");
310 if (reCast<int,IssmDouble>(nodeonSSA[i])){
317 if (!xIsNan<IssmDouble>(spcvx[i])){
321 if (!xIsNan<IssmDouble>(spcvy[i])){
327 else if (reCast<int,IssmDouble>(nodeonFS[i])){
332 if (!xIsNan<IssmDouble>(spcvx[i])){
336 if (!xIsNan<IssmDouble>(spcvy[i])){
340 if (!xIsNan<IssmDouble>(spcvz[i])){
345 else _error_(
"if vertices_type is SSAFS, you shoud have nodeonSSA or nodeonFS");
356 values=xNew<IssmDouble>(Nx);
359 values[j]=spcvx[i*Nx+j];
360 if(!xIsNan<IssmDouble>(values[j]))spcpresent=
true;
367 xDelete<IssmDouble>(values);
380 values=xNew<IssmDouble>(Ny);
383 values[j]=spcvy[i*Ny+j];
384 if(!xIsNan<IssmDouble>(values[j]))spcpresent=
true;
390 xDelete<IssmDouble>(values);
404 values=xNew<IssmDouble>(Nz);
407 values[j]=spcvz[i*Nz+j];
408 if(!xIsNan<IssmDouble>(values[j]))spcpresent=
true;
414 xDelete<IssmDouble>(values);
427 iomodel->
DeleteData(spcvx,
"md.stressbalance.spcvx");
428 iomodel->
DeleteData(spcvy,
"md.stressbalance.spcvy");
429 iomodel->
DeleteData(spcvz,
"md.stressbalance.spcvz");
430 iomodel->
DeleteData(nodeonSSA,
"md.flowequation.borderSSA");
435 iomodel->
DeleteData(vertices_type,
"md.flowequation.vertex_equation");
436 iomodel->
DeleteData(surface,
"md.geometry.surface");
440 xDelete<IssmDouble>(timesx);
441 xDelete<IssmDouble>(timesy);
442 xDelete<IssmDouble>(timesz);
443 xDelete<IssmDouble>(values);
◆ CreateLoads()
void StressbalanceAnalysis::CreateLoads |
( |
Loads * |
loads, |
|
|
IoModel * |
iomodel |
|
) |
| |
|
virtual |
Implements Analysis.
Definition at line 446 of file StressbalanceAnalysis.cpp.
449 const int RIFTINFOSIZE = 12;
453 bool isSSA,isL1L2,isHO,isFS;
454 int numpenalties,numrifts,numriftsegments;
460 iomodel->
FindConstant(&isL1L2,
"md.flowequation.isL1L2");
467 if(!isSSA && !isHO && !isFS && !isL1L2)
return;
473 iomodel->
FetchData(&penalties,&numpenalties,NULL,
"md.stressbalance.vertex_pairing");
475 for(i=0;i<numpenalties;i++){
477 if(iomodel->
my_vertices[reCast<int,IssmDouble>(penalties[2*i+0]-1)]){
480 assert_int=iomodel->
my_vertices[reCast<int,IssmDouble>(penalties[2*i+1]-1)];
_assert_(assert_int);
483 penpair_ids[0]=reCast<int,IssmDouble>(penalties[2*i+0]);
484 penpair_ids[1]=reCast<int,IssmDouble>(penalties[2*i+1]);
493 iomodel->
DeleteData(penalties,
"md.stressbalance.vertex_pairing");
497 iomodel->
FetchData(&riftinfo,&numriftsegments,NULL,
"md.rifts.riftstruct");
498 iomodel->
FetchData(5,
"md.rifts.riftstruct",
"md.geometry.thickness",
"md.geometry.base",
"md.geometry.surface",
"md.mask.ocean_levelset");
499 for(i=0;i<numriftsegments;i++){
500 if(iomodel->
my_elements[reCast<int,IssmDouble>(*(riftinfo+RIFTINFOSIZE*i+2))-1]){
505 iomodel->
DeleteData(5,
"md.rifts.riftstruct",
"md.geometry.thickness",
"md.geometry.base",
"md.geometry.surface",
"md.mask.ocean_levelset");
506 xDelete<IssmDouble>(riftinfo);
◆ CreateNodes()
void StressbalanceAnalysis::CreateNodes |
( |
Nodes * |
nodes, |
|
|
IoModel * |
iomodel, |
|
|
bool |
isamr = false |
|
) |
| |
|
virtual |
Implements Analysis.
Definition at line 509 of file StressbalanceAnalysis.cpp.
512 bool isSSA,isL1L2,isHO,isFS,iscoupling;
513 int finiteelement=-1,approximation=-1;
517 iomodel->
FindConstant(&isL1L2,
"md.flowequation.isL1L2");
522 if(!isSSA & !isL1L2 & !isHO & !isFS)
return;
525 if( (isSSA?1.:0.) + (isL1L2?1.:0.) + (isHO?1.:0.) + (isFS?1.:0.) >1.)
536 iomodel->
FindConstant(&finiteelement,
"md.flowequation.fe_SSA");
544 iomodel->
FindConstant(&finiteelement,
"md.flowequation.fe_HO");
548 iomodel->
FindConstant(&finiteelement,
"md.flowequation.fe_FS");
551 iomodel->
FetchData(3,
"md.flowequation.borderSSA",
"md.flowequation.vertex_equation",
"md.stressbalance.referential");
556 iomodel->
DeleteData(6,
"md.mesh.vertexonbase",
"md.mesh.vertexonsurface",
"md.flowequation.borderSSA",
"md.flowequation.vertex_equation",
557 "md.stressbalance.referential",
"md.flowequation.borderFS");
562 iomodel->
FetchData(6,
"md.mesh.vertexonbase",
"md.mesh.vertexonsurface",
"md.flowequation.borderSSA",
"md.flowequation.borderFS",
563 "md.flowequation.vertex_equation",
"md.stressbalance.referential");
569 approximations[i] = approximation;
574 xDelete<int>(approximations);
576 for(
int i=0;i<nodes->
Size();i++){
578 int sid = node->
Sid();
593 xDelete<int>(approximations);
595 iomodel->
DeleteData(6,
"md.mesh.vertexonbase",
"md.mesh.vertexonsurface",
"md.flowequation.borderSSA",
"md.flowequation.borderFS",
596 "md.flowequation.vertex_equation",
"md.stressbalance.referential");
◆ DofsPerNode()
int StressbalanceAnalysis::DofsPerNode |
( |
int ** |
doflist, |
|
|
int |
domaintype, |
|
|
int |
approximation |
|
) |
| |
|
virtual |
Implements Analysis.
Definition at line 599 of file StressbalanceAnalysis.cpp.
605 switch(approximation){
611 default:
_error_(
"mesh type not supported yet");
619 default:
_error_(
"mesh type not supported yet");
627 default:
_error_(
"mesh type not supported yet");
635 default:
_error_(
"mesh type not supported yet");
640 doftype=xNew<int>(numdofs);
648 doftype=xNew<int>(numdofs);
657 doftype=xNew<int>(numdofs);
◆ UpdateElements()
void StressbalanceAnalysis::UpdateElements |
( |
Elements * |
elements, |
|
|
Inputs2 * |
inputs2, |
|
|
IoModel * |
iomodel, |
|
|
int |
analysis_counter, |
|
|
int |
analysis_type |
|
) |
| |
|
virtual |
Implements Analysis.
Definition at line 672 of file StressbalanceAnalysis.cpp.
675 int materials_type,finiteelement,fe_FS;
676 int approximation,frictionlaw;
677 int FrictionCoupling;
678 int* finiteelement_list=NULL;
679 bool isSSA,isL1L2,isHO,isFS,iscoupling;
680 bool control_analysis;
681 bool dakota_analysis;
686 iomodel->
FindConstant(&isL1L2,
"md.flowequation.isL1L2");
689 iomodel->
FindConstant(&control_analysis,
"md.inversion.iscontrol");
690 iomodel->
FindConstant(&dakota_analysis,
"md.qmu.isdakota");
691 iomodel->
FindConstant(&materials_type,
"md.materials.type");
692 iomodel->
FindConstant(&ismovingfront,
"md.transient.ismovingfront");
696 if(!isSSA & !isL1L2 & !isHO & !isFS)
return;
699 iomodel->
FetchData(1,
"md.flowequation.element_equation");
703 if( (isSSA?1.:0.) + (isL1L2?1.:0.) + (isHO?1.:0.) + (isFS?1.:0.) >1.)
710 if(isSSA) iomodel->
FindConstant(&finiteelement,
"md.flowequation.fe_SSA");
711 else if(isL1L2) finiteelement =
P1Enum;
712 else if(isHO) iomodel->
FindConstant(&finiteelement,
"md.flowequation.fe_HO");
713 else if(isFS) iomodel->
FindConstant(&finiteelement,
"md.flowequation.fe_FS");
715 finiteelement_list[i]=finiteelement;
726 finiteelement_list[i]=
P1Enum;
733 finiteelement_list[i]=finiteelement;
743 element->
Update(inputs2,i,iomodel,analysis_counter,analysis_type,finiteelement_list[i]);
746 if(iomodel->
Data(
"md.flowequation.element_equation")){
765 #ifdef LATERALFRICTION
782 int basalforcing_model;
783 iomodel->
FindConstant(&basalforcing_model,
"md.basalforcings.model");
784 switch(basalforcing_model){
823 iomodel->
FindConstant(&FrictionCoupling,
"md.friction.coupling");
827 if(FrictionCoupling==3){
829 else if(FrictionCoupling==4){
838 iomodel->
FindConstant(&FrictionCoupling,
"md.friction.coupling");
842 if(FrictionCoupling==3){
844 else if(FrictionCoupling==4){
854 iomodel->
FindConstant(&FrictionCoupling,
"md.friction.coupling");
869 iomodel->
FindConstant(&FrictionCoupling,
"md.friction.coupling");
874 if(FrictionCoupling==3){
876 else if(FrictionCoupling==4){
905 _error_(
"friction law "<< frictionlaw <<
" not supported");
908 #ifdef _HAVE_ANDROID_
913 iomodel->
DeleteData(1,
"md.flowequation.element_equation");
914 xDelete<int>(finiteelement_list);
◆ UpdateParameters()
void StressbalanceAnalysis::UpdateParameters |
( |
Parameters * |
parameters, |
|
|
IoModel * |
iomodel, |
|
|
int |
solution_enum, |
|
|
int |
analysis_enum |
|
) |
| |
|
virtual |
Implements Analysis.
Definition at line 916 of file StressbalanceAnalysis.cpp.
921 char** requestedoutputs = NULL;
949 iomodel->
FindConstant(&materials_type,
"md.materials.type");
957 iomodel->
FindConstant(&requestedoutputs,&numoutputs,
"md.stressbalance.requested_outputs");
960 iomodel->
DeleteData(&requestedoutputs,numoutputs,
"md.stressbalance.requested_outputs");
1013 default:
_error_(
"Friction law "<<frictionlaw<<
" not implemented yet");
◆ Core()
void StressbalanceAnalysis::Core |
( |
FemModel * |
femmodel | ) |
|
|
virtual |
Implements Analysis.
Definition at line 1019 of file StressbalanceAnalysis.cpp.
1022 bool isSSA,isL1L2,isHO,isFS;
1023 bool conserve_loads =
true;
1024 int newton,domaintype,fe_FS;
1035 if(isFS && !(isSSA || isHO || isL1L2)){
1040 bool is_schur_cg_solver =
false;
1045 if(solver_type==
FSSolverEnum) is_schur_cg_solver =
true;
1049 if(is_schur_cg_solver)
1060 else if(!isFS && (isSSA || isHO || isL1L2)){
1076 else if ((isSSA || isL1L2 || isHO) && isFS){
◆ CreateDVector()
◆ CreateJacobianMatrix()
◆ CreateKMatrix()
◆ CreatePVector()
◆ GetSolutionFromInputs()
◆ GetSolutionFromInputsHoriz()
Definition at line 1189 of file StressbalanceAnalysis.cpp.
1192 int domaintype,dim,approximation,dofpernode;
1193 int* doflist = NULL;
1206 int numdof = numnodes*dofpernode;
1211 IssmDouble* values = xNew<IssmDouble>(numdof);
1220 for(
int i=0;i<numnodes;i++){
1225 values[i*dofpernode+0]=vx;
1228 values[i*dofpernode+1]=vy;
1236 xDelete<IssmDouble>(values);
1237 xDelete<int>(doflist);
◆ GradientJ()
void StressbalanceAnalysis::GradientJ |
( |
Vector< IssmDouble > * |
gradient, |
|
|
Element * |
element, |
|
|
int |
control_type, |
|
|
int |
control_index |
|
) |
| |
|
virtual |
◆ InputUpdateFromSolution()
void StressbalanceAnalysis::InputUpdateFromSolution |
( |
IssmDouble * |
solution, |
|
|
Element * |
element |
|
) |
| |
|
virtual |
◆ UpdateConstraints()
void StressbalanceAnalysis::UpdateConstraints |
( |
FemModel * |
femmodel | ) |
|
|
virtual |
◆ CreateJacobianMatrixSSA()
Definition at line 1279 of file StressbalanceAnalysis.cpp.
1289 basalelement = element;
1292 if(!element->
IsOnBase())
return NULL;
1312 IssmDouble* dbasis = xNew<IssmDouble>(2*numnodes);
1322 for(
int ig=gauss->
begin();ig<gauss->end();ig++){
1329 basalelement->
StrainRateSSA(&epsilon[0],xyz_list,gauss,vx_input,vy_input);
1331 eps1[0]=2*epsilon[0]+epsilon[1]; eps2[0]=epsilon[2];
1332 eps1[1]=epsilon[2]; eps2[1]=epsilon[0]+2*epsilon[1];
1334 for(
int i=0;i<numnodes;i++){
1335 for(
int j=0;j<numnodes;j++){
1336 eps1dotdphii=eps1[0]*dbasis[0*numnodes+i]+eps1[1]*dbasis[1*numnodes+i];
1337 eps1dotdphij=eps1[0]*dbasis[0*numnodes+j]+eps1[1]*dbasis[1*numnodes+j];
1338 eps2dotdphii=eps2[0]*dbasis[0*numnodes+i]+eps2[1]*dbasis[1*numnodes+i];
1339 eps2dotdphij=eps2[0]*dbasis[0*numnodes+j]+eps2[1]*dbasis[1*numnodes+j];
1341 Ke->
values[2*numnodes*(2*i+0)+2*j+0]+=gauss->
weight*Jdet*2.*mu_prime*thickness*eps1dotdphij*eps1dotdphii;
1342 Ke->
values[2*numnodes*(2*i+0)+2*j+1]+=gauss->
weight*Jdet*2.*mu_prime*thickness*eps2dotdphij*eps1dotdphii;
1343 Ke->
values[2*numnodes*(2*i+1)+2*j+0]+=gauss->
weight*Jdet*2.*mu_prime*thickness*eps1dotdphij*eps2dotdphii;
1344 Ke->
values[2*numnodes*(2*i+1)+2*j+1]+=gauss->
weight*Jdet*2.*mu_prime*thickness*eps2dotdphij*eps2dotdphii;
1354 xDelete<IssmDouble>(xyz_list);
1355 xDelete<IssmDouble>(dbasis);
◆ CreateKMatrixSSA()
◆ CreateKMatrixSSAFriction()
Definition at line 1401 of file StressbalanceAnalysis.cpp.
1408 bool mainlyfloating;
1409 int friction_style,point1;
1413 Gauss* gauss = NULL;
1426 int numdof = numnodes*dim;
1430 IssmDouble* basis = xNew<IssmDouble>(numnodes);
1436 Input2* gllevelset_input = NULL;
1445 element->
GetGroundedPart(&point1,&fraction1,&fraction2,&mainlyfloating);
1446 gauss = element->
NewGauss(point1,fraction1,fraction2,mainlyfloating,2);
1453 for(
int ig=gauss->
begin();ig<gauss->end();ig++){
1460 if(gllevelset<0.) alpha2=0.;
1463 if (phi<0.99999999) alpha2=0.;
1465 else _error_(
"friction interpolation "<<
EnumToStringx(friction_style)<<
" not implemented yet");
1471 for(
int i=0;i<numnodes;i++){
1472 for(
int j=0;j<numnodes;j++){
1473 Ke->
values[2*i*2*numnodes+2*j] += alpha2*gauss->
weight*Jdet*basis[i]*basis[j];
1474 Ke->
values[(2*i+1)*2*numnodes+2*j+1] += alpha2*gauss->
weight*Jdet*basis[i]*basis[j];
1479 for(
int i=0;i<numnodes;i++){
1480 for(
int j=0;j<numnodes;j++){
1481 Ke->
values[i*numnodes+j] += alpha2*gauss->
weight*Jdet*basis[i]*basis[j];
1493 xDelete<IssmDouble>(xyz_list);
1494 xDelete<IssmDouble>(basis);
◆ CreateKMatrixSSALateralFriction()
Definition at line 1497 of file StressbalanceAnalysis.cpp.
1520 int numdof = numnodes*dim;
1524 IssmDouble* B = xNew<IssmDouble>(dim*numdof);
1525 IssmDouble* D = xNewZeroInit<IssmDouble>(dim*dim);
1533 Gauss* gauss=element->
NewGauss(xyz_list,xyz_list_boundary,3);
1534 for(
int ig=gauss->
begin();ig<gauss->end();ig++){
1544 for(
int i=0;i<dim;i++) D[i*dim+i]=alpha2*gauss->
weight*Jdet;
1557 xDelete<IssmDouble>(xyz_list);
1558 xDelete<IssmDouble>(xyz_list_boundary);
1559 xDelete<IssmDouble>(B);
1560 xDelete<IssmDouble>(D);
◆ CreateKMatrixSSAViscous()
Definition at line 1563 of file StressbalanceAnalysis.cpp.
1584 int numdof = numnodes*dim;
1588 IssmDouble* dbasis = xNew<IssmDouble>(2*numnodes);
1601 for(
int ig=gauss->
begin();ig<gauss->end();ig++){
1611 for(
int i=0;i<numnodes;i++){
1612 for(
int j=0;j<numnodes;j++){
1613 Ke->
values[2*i*2*numnodes+2*j] += gauss->
weight*Jdet*viscosity*thickness*(
1614 4.*dbasis[0*numnodes+j]*dbasis[0*numnodes+i] + dbasis[1*numnodes+j]*dbasis[1*numnodes+i]
1616 Ke->
values[2*i*2*numnodes+2*j+1] += gauss->
weight*Jdet*viscosity*thickness*(
1617 2.*dbasis[1*numnodes+j]*dbasis[0*numnodes+i] + dbasis[0*numnodes+j]*dbasis[1*numnodes+i]
1619 Ke->
values[(2*i+1)*2*numnodes+2*j] += gauss->
weight*Jdet*viscosity*thickness*(
1620 2.*dbasis[0*numnodes+j]*dbasis[1*numnodes+i] + dbasis[1*numnodes+j]*dbasis[0*numnodes+i]
1622 Ke->
values[(2*i+1)*2*numnodes+2*j+1] += gauss->
weight*Jdet*viscosity*thickness*(
1623 dbasis[0*numnodes+j]*dbasis[0*numnodes+i] + 4.*dbasis[1*numnodes+j]*dbasis[1*numnodes+i]
1629 for(
int i=0;i<numnodes;i++){
1630 for(
int j=0;j<numnodes;j++){
1631 Ke->
values[i*numnodes+j] += gauss->
weight*Jdet*viscosity*thickness*(
1632 4.*dbasis[0*numnodes+j]*dbasis[0*numnodes+i]
1644 xDelete<IssmDouble>(xyz_list);
1645 xDelete<IssmDouble>(dbasis);
◆ CreatePVectorSSA()
◆ CreatePVectorSSAFront()
Definition at line 1740 of file StressbalanceAnalysis.cpp.
1750 IssmDouble Jdet,thickness,base,sealevel,water_pressure,ice_pressure;
1751 IssmDouble surface_under_water,base_under_water,pressure;
1770 IssmDouble* basis = xNew<IssmDouble>(numnodes);
1785 for(
int ig=gauss->
begin();ig<gauss->end();ig++){
1794 surface_under_water =
min(0.,thickness+base-sealevel);
1795 base_under_water =
min(0.,base-sealevel);
1796 water_pressure = 1.0/2.0*gravity*rho_water*(surface_under_water*surface_under_water - base_under_water*base_under_water);
1797 ice_pressure = 1.0/2.0*gravity*rho_ice*thickness*thickness;
1798 pressure = ice_pressure + water_pressure;
1800 for(
int i=0;i<numnodes;i++){
1801 pe->
values[dim*i+0]+= pressure*Jdet*gauss->
weight*normal[0]*basis[i];
1802 if(dim==2) pe->
values[dim*i+1]+= pressure*Jdet*gauss->
weight*normal[1]*basis[i];
1810 xDelete<IssmDouble>(xyz_list);
1811 xDelete<IssmDouble>(xyz_list_front);
1812 xDelete<IssmDouble>(basis);
◆ CreatePVectorSSADrivingStress()
Definition at line 1681 of file StressbalanceAnalysis.cpp.
1705 IssmDouble* basis = xNew<IssmDouble>(numnodes);
1715 for(
int ig=gauss->
begin();ig<gauss->end();ig++){
1725 for(
int i=0;i<numnodes;i++){
1726 pe->
values[i*dim+0]+=-rhog*thickness*slope[0]*Jdet*gauss->
weight*basis[i];
1727 if(dim==2) pe->
values[i*dim+1]+=-rhog*thickness*slope[1]*Jdet*gauss->
weight*basis[i];
1735 xDelete<IssmDouble>(xyz_list);
1736 xDelete<IssmDouble>(basis);
◆ GetBSSA()
Definition at line 1816 of file StressbalanceAnalysis.cpp.
1833 IssmDouble* dbasis=xNew<IssmDouble>(dim*numnodes);
1838 for(
int i=0;i<numnodes;i++){
1839 B[2*numnodes*0+2*i+0] = dbasis[0*numnodes+i];
1840 B[2*numnodes*0+2*i+1] = 0.;
1841 B[2*numnodes*1+2*i+0] = 0.;
1842 B[2*numnodes*1+2*i+1] = dbasis[1*numnodes+i];
1843 B[2*numnodes*2+2*i+0] = .5*dbasis[1*numnodes+i];
1844 B[2*numnodes*2+2*i+1] = .5*dbasis[0*numnodes+i];
1848 for(
int i=0;i<numnodes;i++){
1854 xDelete<IssmDouble>(dbasis);
◆ GetBSSAFriction()
Definition at line 1856 of file StressbalanceAnalysis.cpp.
1872 IssmDouble* basis=xNew<IssmDouble>(numnodes);
1877 for(
int i=0;i<numnodes;i++){
1878 B[2*numnodes*0+2*i+0] = basis[i];
1879 B[2*numnodes*0+2*i+1] = 0.;
1880 B[2*numnodes*1+2*i+0] = 0.;
1881 B[2*numnodes*1+2*i+1] = basis[i];
1885 for(
int i=0;i<numnodes;i++){
1891 xDelete<IssmDouble>(basis);
◆ GetBSSAprime()
Definition at line 1893 of file StressbalanceAnalysis.cpp.
1910 IssmDouble* dbasis=xNew<IssmDouble>(dim*numnodes);
1915 for(
int i=0;i<numnodes;i++){
1916 Bprime[2*numnodes*0+2*i+0] = 2.*dbasis[0*numnodes+i];
1917 Bprime[2*numnodes*0+2*i+1] = dbasis[1*numnodes+i];
1918 Bprime[2*numnodes*1+2*i+0] = dbasis[0*numnodes+i];
1919 Bprime[2*numnodes*1+2*i+1] = 2.*dbasis[1*numnodes+i];
1920 Bprime[2*numnodes*2+2*i+0] = dbasis[1*numnodes+i];
1921 Bprime[2*numnodes*2+2*i+1] = dbasis[0*numnodes+i];
1925 for(
int i=0;i<numnodes;i++){
1926 Bprime[i] = 2.*dbasis[i];
1931 xDelete<IssmDouble>(dbasis);
◆ InputUpdateFromSolutionSSA()
void StressbalanceAnalysis::InputUpdateFromSolutionSSA |
( |
IssmDouble * |
solution, |
|
|
Element * |
element |
|
) |
| |
Definition at line 1933 of file StressbalanceAnalysis.cpp.
1935 int i,dim,domaintype;
1943 IssmDouble* pressure = xNew<IssmDouble>(numvertices);
1944 IssmDouble* thickness = xNew<IssmDouble>(numvertices);
1945 IssmDouble* surface = xNew<IssmDouble>(numvertices);
1953 for(i=0;i<numvertices;i++) pressure[i]=rho_ice*g*thickness[i];
1959 for(i=0;i<numvertices;i++) pressure[i]=rho_ice*g*(surface[i]-xyz_list[i*3+2]);
1965 for(i=0;i<numvertices;i++) pressure[i]=rho_ice*g*(surface[i]-xyz_list[i*3+1]);
1971 xDelete<IssmDouble>(pressure);
1972 xDelete<IssmDouble>(thickness);
1973 xDelete<IssmDouble>(surface);
1978 basalelement = element;
1981 if(!element->
IsOnBase()){xDelete<IssmDouble>(xyz_list);
return;}
1989 int numdof = numnodes*dim;
1993 IssmDouble* values = xNew<IssmDouble>(numdof);
1997 IssmDouble* vel = xNew<IssmDouble>(numnodes);
2000 for(i=0;i<numdof;i++) values[i]=solution[doflist[i]];
2006 for(i=0;i<numnodes;i++){
2007 vx[i]=values[i*dim+0];
2008 if(xIsNan<IssmDouble>(vx[i]))
_error_(
"NaN found in solution vector");
2009 if(xIsInf<IssmDouble>(vx[i]))
_error_(
"Inf found in solution vector");
2012 vy[i]=values[i*dim+1];
2013 if(xIsNan<IssmDouble>(vy[i]))
_error_(
"NaN found in solution vector");
2014 if(xIsInf<IssmDouble>(vy[i]))
_error_(
"Inf found in solution vector");
2021 for(i=0;i<numnodes;i++) vel[i]=sqrt(vx[i]*vx[i] + vy[i]*vy[i] + vz[i]*vz[i]);
2025 for(i=0;i<numnodes;i++) vel[i]=sqrt(vx[i]*vx[i] + vy[i]*vy[i]);
2034 xDelete<IssmDouble>(vel);
2035 xDelete<IssmDouble>(vz);
2036 xDelete<IssmDouble>(vy);
2037 xDelete<IssmDouble>(vx);
2038 xDelete<IssmDouble>(values);
2039 xDelete<IssmDouble>(xyz_list);
2040 xDelete<int>(doflist);
◆ CreateKMatrixL1L2()
◆ CreateKMatrixL1L2Friction()
◆ CreateKMatrixL1L2Viscous()
Definition at line 2070 of file StressbalanceAnalysis.cpp.
2081 int numdof = numnodes*2;
2086 IssmDouble* Bprime = xNew<IssmDouble>(3*numdof);
2087 IssmDouble* D = xNewZeroInit<IssmDouble>(3*3);
2098 for(
int ig=gauss->
begin();ig<gauss->end();ig++){
2104 this->
GetBSSA(B,basalelement,2,xyz_list,gauss_base);
2105 this->
GetBSSAprime(Bprime,basalelement,2,xyz_list,gauss_base);
2109 for(
int i=0;i<3;i++) D[i*3+i]=2*viscosity*gauss->
weight*Jdet;
2124 xDelete<IssmDouble>(xyz_list);
2125 xDelete<IssmDouble>(D);
2126 xDelete<IssmDouble>(Bprime);
2127 xDelete<IssmDouble>(B);
◆ CreatePVectorL1L2()
◆ CreatePVectorL1L2Front()
Definition at line 2209 of file StressbalanceAnalysis.cpp.
2215 IssmDouble Jdet,thickness,bed,sealevel,water_pressure,ice_pressure;
2216 IssmDouble surface_under_water,base_under_water,pressure;
2226 IssmDouble* basis = xNew<IssmDouble>(numnodes);
2241 for(
int ig=gauss->
begin();ig<gauss->end();ig++){
2250 surface_under_water =
min(0.,thickness+bed-sealevel);
2251 base_under_water =
min(0.,bed-sealevel);
2252 water_pressure = 1.0/2.0*gravity*rho_water*(surface_under_water*surface_under_water - base_under_water*base_under_water);
2253 ice_pressure = 1.0/2.0*gravity*rho_ice*thickness*thickness;
2254 pressure = ice_pressure + water_pressure;
2256 for (
int i=0;i<numnodes;i++){
2257 pe->
values[2*i+0]+= pressure*Jdet*gauss->
weight*normal[0]*basis[i];
2258 pe->
values[2*i+1]+= pressure*Jdet*gauss->
weight*normal[1]*basis[i];
2266 xDelete<IssmDouble>(xyz_list);
2267 xDelete<IssmDouble>(xyz_list_front);
2268 xDelete<IssmDouble>(basis);
◆ CreatePVectorL1L2DrivingStress()
Definition at line 2163 of file StressbalanceAnalysis.cpp.
2174 IssmDouble* basis = xNew<IssmDouble>(numnodes);
2184 for(
int ig=gauss->
begin();ig<gauss->end();ig++){
2194 for(
int i=0;i<numnodes;i++){
2195 pe->
values[i*2+0]+=-rhog*thickness*slope[0]*Jdet*gauss->
weight*basis[i];
2196 pe->
values[i*2+1]+=-rhog*thickness*slope[1]*Jdet*gauss->
weight*basis[i];
2204 xDelete<IssmDouble>(xyz_list);
2205 xDelete<IssmDouble>(basis);
◆ InputUpdateFromSolutionL1L2()
void StressbalanceAnalysis::InputUpdateFromSolutionL1L2 |
( |
IssmDouble * |
solution, |
|
|
Element * |
element |
|
) |
| |
Definition at line 2272 of file StressbalanceAnalysis.cpp.
2274 int i,dim,domaintype;
2282 IssmDouble* pressure = xNew<IssmDouble>(numvertices);
2283 IssmDouble* thickness = xNew<IssmDouble>(numvertices);
2284 IssmDouble* surface = xNew<IssmDouble>(numvertices);
2292 for(i=0;i<numvertices;i++) pressure[i]=rho_ice*g*thickness[i];
2297 for(i=0;i<numvertices;i++) pressure[i]=rho_ice*g*(surface[i]-xyz_list[i*3+2]);
2300 xDelete<IssmDouble>(pressure);
2301 xDelete<IssmDouble>(thickness);
2302 xDelete<IssmDouble>(surface);
2307 basalelement = element;
2310 if(!element->
IsOnBase()){xDelete<IssmDouble>(xyz_list);
return;}
2318 int numdof = numnodes*2;
2322 IssmDouble* values = xNew<IssmDouble>(numdof);
2326 IssmDouble* vel = xNew<IssmDouble>(numnodes);
2329 for(i=0;i<numdof;i++) values[i]=solution[doflist[i]];
2336 for(i=0;i<numnodes;i++){
2337 vx[i]=values[i*2+0];
2338 vy[i]=values[i*2+1];
2341 if(xIsNan<IssmDouble>(vx[i]))
_error_(
"NaN found in solution vector");
2342 if(xIsInf<IssmDouble>(vx[i]))
_error_(
"Inf found in solution vector");
2343 if(xIsNan<IssmDouble>(vy[i]))
_error_(
"NaN found in solution vector");
2344 if(xIsInf<IssmDouble>(vy[i]))
_error_(
"Inf found in solution vector");
2349 for(i=0;i<numnodes;i++) vel[i]=sqrt(vx[i]*vx[i] + vy[i]*vy[i] + vz[i]*vz[i]);
2357 xDelete<IssmDouble>(vel);
2358 xDelete<IssmDouble>(vz);
2359 xDelete<IssmDouble>(vy);
2360 xDelete<IssmDouble>(vx);
2361 xDelete<IssmDouble>(values);
2362 xDelete<IssmDouble>(xyz_list);
2363 xDelete<int>(doflist);
◆ CreateJacobianMatrixHO()
Definition at line 2368 of file StressbalanceAnalysis.cpp.
2384 IssmDouble* dbasis = xNew<IssmDouble>(3*numnodes);
2393 for(
int ig=gauss->
begin();ig<gauss->end();ig++){
2399 element->
StrainRateHO(&epsilon[0],xyz_list,gauss,vx_input,vy_input);
2401 eps1[0]=2*epsilon[0]+epsilon[1]; eps2[0]=epsilon[2];
2402 eps1[1]=epsilon[2]; eps2[1]=epsilon[0]+2*epsilon[1];
2403 eps1[2]=epsilon[3]; eps2[2]=epsilon[4];
2405 for(
int i=0;i<numnodes;i++){
2406 for(
int j=0;j<numnodes;j++){
2407 eps1dotdphii=eps1[0]*dbasis[0*numnodes+i]+eps1[1]*dbasis[1*numnodes+i]+eps1[2]*dbasis[2*numnodes+i];
2408 eps1dotdphij=eps1[0]*dbasis[0*numnodes+j]+eps1[1]*dbasis[1*numnodes+j]+eps1[2]*dbasis[2*numnodes+j];
2409 eps2dotdphii=eps2[0]*dbasis[0*numnodes+i]+eps2[1]*dbasis[1*numnodes+i]+eps2[2]*dbasis[2*numnodes+i];
2410 eps2dotdphij=eps2[0]*dbasis[0*numnodes+j]+eps2[1]*dbasis[1*numnodes+j]+eps2[2]*dbasis[2*numnodes+j];
2412 Ke->
values[2*numnodes*(2*i+0)+2*j+0]+=gauss->
weight*Jdet*2.*mu_prime*eps1dotdphij*eps1dotdphii;
2413 Ke->
values[2*numnodes*(2*i+0)+2*j+1]+=gauss->
weight*Jdet*2.*mu_prime*eps2dotdphij*eps1dotdphii;
2414 Ke->
values[2*numnodes*(2*i+1)+2*j+0]+=gauss->
weight*Jdet*2.*mu_prime*eps1dotdphij*eps2dotdphii;
2415 Ke->
values[2*numnodes*(2*i+1)+2*j+1]+=gauss->
weight*Jdet*2.*mu_prime*eps2dotdphij*eps2dotdphii;
2425 xDelete<IssmDouble>(xyz_list);
2426 xDelete<IssmDouble>(dbasis);
◆ CreateKMatrixHO()
◆ CreateKMatrixHOFriction()
Definition at line 2444 of file StressbalanceAnalysis.cpp.
2453 bool mainlyfloating;
2454 int friction_style,point1;
2458 Gauss* gauss = NULL;
2465 int numdof = numnodes*(dim-1);
2469 IssmDouble* basis = xNew<IssmDouble>(numnodes);
2474 Input2* gllevelset_input = NULL;
2483 element->
GetGroundedPart(&point1,&fraction1,&fraction2,&mainlyfloating);
2484 gauss = element->
NewGauss(point1,fraction1,fraction2,mainlyfloating,2);
2491 for(
int ig=gauss->
begin();ig<gauss->end();ig++){
2498 if(gllevelset<0.) alpha2=0.;
2501 if (phi<0.99999999) alpha2=0.;
2503 else _error_(
"friction interpolation "<<
EnumToStringx(friction_style)<<
" not implemented yet");
2509 for(
int i=0;i<numnodes;i++){
2510 for(
int j=0;j<numnodes;j++){
2511 Ke->values[2*i*2*numnodes+2*j] += alpha2*gauss->
weight*Jdet*basis[i]*basis[j];
2512 Ke->values[(2*i+1)*2*numnodes+2*j+1] += alpha2*gauss->
weight*Jdet*basis[i]*basis[j];
2517 for(
int i=0;i<numnodes;i++){
2518 for(
int j=0;j<numnodes;j++){
2519 Ke->values[i*numnodes+j] += alpha2*gauss->
weight*Jdet*basis[i]*basis[j];
2531 xDelete<IssmDouble>(xyz_list_base);
2532 xDelete<IssmDouble>(basis);
◆ CreateKMatrixHOViscous()
Definition at line 2535 of file StressbalanceAnalysis.cpp.
2550 int numdof = numnodes*(dim-1);
2554 IssmDouble* dbasis = xNew<IssmDouble>(dim*numnodes);
2566 for(
int ig=gauss->
begin();ig<gauss->end();ig++){
2574 for(
int i=0;i<numnodes;i++){
2575 for(
int j=0;j<numnodes;j++){
2576 Ke->values[2*i*2*numnodes+2*j] += gauss->
weight*Jdet*viscosity*(
2577 4.*dbasis[0*numnodes+j]*dbasis[0*numnodes+i] + dbasis[1*numnodes+j]*dbasis[1*numnodes+i] + dbasis[2*numnodes+j]*dbasis[2*numnodes+i]
2579 Ke->values[2*i*2*numnodes+2*j+1] += gauss->
weight*Jdet*viscosity*(
2580 2.*dbasis[1*numnodes+j]*dbasis[0*numnodes+i] + dbasis[0*numnodes+j]*dbasis[1*numnodes+i]
2582 Ke->values[(2*i+1)*2*numnodes+2*j] += gauss->
weight*Jdet*viscosity*(
2583 2.*dbasis[0*numnodes+j]*dbasis[1*numnodes+i] + dbasis[1*numnodes+j]*dbasis[0*numnodes+i]
2585 Ke->values[(2*i+1)*2*numnodes+2*j+1] += gauss->
weight*Jdet*viscosity*(
2586 dbasis[0*numnodes+j]*dbasis[0*numnodes+i] + 4.*dbasis[1*numnodes+j]*dbasis[1*numnodes+i] + dbasis[2*numnodes+j]*dbasis[2*numnodes+i]
2592 for(
int i=0;i<numnodes;i++){
2593 for(
int j=0;j<numnodes;j++){
2594 Ke->values[i*numnodes+j] += gauss->
weight*Jdet*viscosity*(
2595 4.*dbasis[0*numnodes+j]*dbasis[0*numnodes+i] + dbasis[1*numnodes+j]*dbasis[1*numnodes+i]
2607 xDelete<IssmDouble>(xyz_list);
2608 xDelete<IssmDouble>(dbasis);
◆ CreatePVectorHO()
◆ CreatePVectorHOFront()
Definition at line 2733 of file StressbalanceAnalysis.cpp.
2743 IssmDouble Jdet,surface,sealevel,z,water_pressure,ice_pressure;
2744 IssmDouble surface_under_water,base_under_water,pressure;
2748 Gauss* gauss = NULL;
2759 IssmDouble* basis = xNew<IssmDouble>(numnodes);
2772 IssmDouble zmax=xyz_list[0*3+(dim-1)];
for(
int i=1;i<numvertices;i++)
if(xyz_list[i*3+(dim-1)]>zmax) zmax=xyz_list[i*3+(dim-1)];
2773 IssmDouble zmin=xyz_list[0*3+(dim-1)];
for(
int i=1;i<numvertices;i++)
if(xyz_list[i*3+(dim-1)]<zmin) zmin=xyz_list[i*3+(dim-1)];
2774 if(zmax>0. && zmin<0.) gauss=element->
NewGauss(xyz_list,xyz_list_front,3,10);
2775 else gauss=element->
NewGauss(xyz_list,xyz_list_front,3,3);
2778 for(
int ig=gauss->
begin();ig<gauss->
end();ig++){
2783 if(dim==3) z=element->
GetZcoord(xyz_list,gauss);
2784 else z=element->
GetYcoord(xyz_list,gauss);
2788 water_pressure = rho_water*gravity*
min(0.,z-sealevel);
2789 ice_pressure = rho_ice*gravity*(surface-z);
2790 pressure = ice_pressure + water_pressure;
2792 for (
int i=0;i<numnodes;i++){
2793 pe->
values[(dim-1)*i+0]+= pressure*Jdet*gauss->
weight*normal[0]*basis[i];
2794 if(dim==3) pe->
values[(dim-1)*i+1]+= pressure*Jdet*gauss->
weight*normal[1]*basis[i];
2802 xDelete<IssmDouble>(basis);
2803 xDelete<IssmDouble>(xyz_list);
2804 xDelete<IssmDouble>(xyz_list_front);
◆ CreatePVectorHODrivingStress()
Definition at line 2684 of file StressbalanceAnalysis.cpp.
2702 IssmDouble* basis = xNew<IssmDouble>(numnodes);
2711 for(
int ig=gauss->
begin();ig<gauss->end();ig++){
2718 for(
int i=0;i<numnodes;i++){
2719 pe->
values[i*(dim-1)+0]+=-rhog*slope[0]*Jdet*gauss->
weight*basis[i];
2720 if(dim==3) pe->
values[i*(dim-1)+1]+=-rhog*slope[1]*Jdet*gauss->
weight*basis[i];
2728 xDelete<IssmDouble>(basis);
2729 xDelete<IssmDouble>(xyz_list);
◆ GetBHO()
Definition at line 2808 of file StressbalanceAnalysis.cpp.
2828 IssmDouble* dbasis=xNew<IssmDouble>(dim*numnodes);
2833 for(
int i=0;i<numnodes;i++){
2834 B[numnodes*0+i] = dbasis[0*numnodes+i];
2835 B[numnodes*1+i] = .5*dbasis[1*numnodes+i];
2839 for(
int i=0;i<numnodes;i++){
2840 B[2*numnodes*0+2*i+0] = dbasis[0*numnodes+i];
2841 B[2*numnodes*0+2*i+1] = 0.;
2842 B[2*numnodes*1+2*i+0] = 0.;
2843 B[2*numnodes*1+2*i+1] = dbasis[1*numnodes+i];
2844 B[2*numnodes*2+2*i+0] = .5*dbasis[1*numnodes+i];
2845 B[2*numnodes*2+2*i+1] = .5*dbasis[0*numnodes+i];
2846 B[2*numnodes*3+2*i+0] = .5*dbasis[2*numnodes+i];
2847 B[2*numnodes*3+2*i+1] = 0.;
2848 B[2*numnodes*4+2*i+0] = 0.;
2849 B[2*numnodes*4+2*i+1] = .5*dbasis[2*numnodes+i];
2854 xDelete<IssmDouble>(dbasis);
◆ GetBHOFriction()
Definition at line 2856 of file StressbalanceAnalysis.cpp.
2872 IssmDouble* basis=xNew<IssmDouble>(numnodes);
2877 for(
int i=0;i<numnodes;i++){
2878 B[2*numnodes*0+2*i+0] = basis[i];
2879 B[2*numnodes*0+2*i+1] = 0.;
2880 B[2*numnodes*1+2*i+0] = 0.;
2881 B[2*numnodes*1+2*i+1] = basis[i];
2885 for(
int i=0;i<numnodes;i++){
2891 xDelete<IssmDouble>(basis);
◆ GetBHOprime()
Definition at line 2893 of file StressbalanceAnalysis.cpp.
2910 IssmDouble* dbasis=xNew<IssmDouble>(dim*numnodes);
2915 for(
int i=0;i<numnodes;i++){
2916 Bprime[2*numnodes*0+2*i+0] = 2.*dbasis[0*numnodes+i];
2917 Bprime[2*numnodes*0+2*i+1] = dbasis[1*numnodes+i];
2918 Bprime[2*numnodes*1+2*i+0] = dbasis[0*numnodes+i];
2919 Bprime[2*numnodes*1+2*i+1] = 2.*dbasis[1*numnodes+i];
2920 Bprime[2*numnodes*2+2*i+0] = dbasis[1*numnodes+i];
2921 Bprime[2*numnodes*2+2*i+1] = dbasis[0*numnodes+i];
2922 Bprime[2*numnodes*3+2*i+0] = dbasis[2*numnodes+i];
2923 Bprime[2*numnodes*3+2*i+1] = 0.;
2924 Bprime[2*numnodes*4+2*i+0] = 0.;
2925 Bprime[2*numnodes*4+2*i+1] = dbasis[2*numnodes+i];
2929 for(
int i=0;i<numnodes;i++){
2930 Bprime[numnodes*0+i] = 2.*dbasis[0*numnodes+i];
2931 Bprime[numnodes*1+i] = dbasis[1*numnodes+i];
2936 xDelete<IssmDouble>(dbasis);
◆ InputUpdateFromSolutionHO()
void StressbalanceAnalysis::InputUpdateFromSolutionHO |
( |
IssmDouble * |
solution, |
|
|
Element * |
element |
|
) |
| |
Definition at line 2938 of file StressbalanceAnalysis.cpp.
2940 int i,domaintype,dim;
2946 IssmDouble* pressure = xNew<IssmDouble>(numvertices);
2947 IssmDouble* surface = xNew<IssmDouble>(numvertices);
2956 for(i=0;i<numvertices;i++) pressure[i]=rho_ice*g*(surface[i]-xyz_list[i*3+2]);
2962 for(i=0;i<numvertices;i++) pressure[i]=rho_ice*g*(surface[i]-xyz_list[i*3+1]);
2968 xDelete<IssmDouble>(pressure);
2969 xDelete<IssmDouble>(surface);
2973 int numdof = numnodes*(dim-1);
2977 IssmDouble* values = xNew<IssmDouble>(numdof);
2981 IssmDouble* vel = xNew<IssmDouble>(numnodes);
2984 for(i=0;i<numdof;i++) values[i]=solution[doflist[i]];
2990 for(i=0;i<numnodes;i++){
2991 vx[i]=values[i*(dim-1)+0];
2992 if(xIsNan<IssmDouble>(vx[i]))
_error_(
"NaN found in solution vector");
2993 if(xIsInf<IssmDouble>(vx[i]))
_error_(
"Inf found in solution vector");
2995 vy[i]=values[i*(dim-1)+1];
2996 if(xIsNan<IssmDouble>(vy[i]))
_error_(
"NaN found in solution vector");
2997 if(xIsInf<IssmDouble>(vy[i]))
_error_(
"Inf found in solution vector");
3004 for(i=0;i<numnodes;i++) vel[i]=sqrt(vx[i]*vx[i] + vy[i]*vy[i] + vz[i]*vz[i]);
3008 for(i=0;i<numnodes;i++) vel[i]=sqrt(vx[i]*vx[i] + vy[i]*vy[i]);
3017 xDelete<IssmDouble>(vel);
3018 xDelete<IssmDouble>(vz);
3019 xDelete<IssmDouble>(vy);
3020 xDelete<IssmDouble>(vx);
3021 xDelete<IssmDouble>(values);
3022 xDelete<IssmDouble>(xyz_list);
3023 xDelete<int>(doflist);
◆ CreateDVectorFS()
◆ CreateJacobianMatrixFS()
Definition at line 3051 of file StressbalanceAnalysis.cpp.
3067 int numdof = vnumnodes*3 + pnumnodes*1;
3070 int* cs_list = xNew<int>(vnumnodes+pnumnodes);
3071 for(i=0;i<vnumnodes;i++) cs_list[i] =
XYZEnum;
3072 for(i=0;i<pnumnodes;i++) cs_list[vnumnodes+i] =
PressureEnum;
3076 IssmDouble* dbasis = xNew<IssmDouble>(3*vnumnodes);
3086 for(
int ig=gauss->
begin();ig<gauss->end();ig++){
3096 element->
StrainRateHO(&epsilon[0],xyz_list,gauss,vx_input,vy_input);
3097 eps1[0]=epsilon[0]; eps2[0]=epsilon[2]; eps3[0]=epsilon[3];
3098 eps1[1]=epsilon[2]; eps2[1]=epsilon[1]; eps3[1]=epsilon[4];
3099 eps1[2]=epsilon[3]; eps2[2]=epsilon[4]; eps3[2]= -epsilon[0] -epsilon[1];
3102 for(i=0;i<vnumnodes;i++){
3103 for(j=0;j<vnumnodes;j++){
3104 eps1dotdphii=eps1[0]*dbasis[0*vnumnodes+i]+eps1[1]*dbasis[1*vnumnodes+i]+eps1[2]*dbasis[2*vnumnodes+i];
3105 eps1dotdphij=eps1[0]*dbasis[0*vnumnodes+j]+eps1[1]*dbasis[1*vnumnodes+j]+eps1[2]*dbasis[2*vnumnodes+j];
3106 eps2dotdphii=eps2[0]*dbasis[0*vnumnodes+i]+eps2[1]*dbasis[1*vnumnodes+i]+eps2[2]*dbasis[2*vnumnodes+i];
3107 eps2dotdphij=eps2[0]*dbasis[0*vnumnodes+j]+eps2[1]*dbasis[1*vnumnodes+j]+eps2[2]*dbasis[2*vnumnodes+j];
3108 eps3dotdphii=eps3[0]*dbasis[0*vnumnodes+i]+eps3[1]*dbasis[1*vnumnodes+i]+eps3[2]*dbasis[2*vnumnodes+i];
3109 eps3dotdphij=eps3[0]*dbasis[0*vnumnodes+j]+eps3[1]*dbasis[1*vnumnodes+j]+eps3[2]*dbasis[2*vnumnodes+j];
3111 Ke->
values[numdof*(3*i+0)+3*j+0]+=gauss->
weight*Jdet*2*mu_prime*eps1dotdphij*eps1dotdphii;
3112 Ke->
values[numdof*(3*i+0)+3*j+1]+=gauss->
weight*Jdet*2*mu_prime*eps2dotdphij*eps1dotdphii;
3113 Ke->
values[numdof*(3*i+0)+3*j+2]+=gauss->
weight*Jdet*2*mu_prime*eps3dotdphij*eps1dotdphii;
3115 Ke->
values[numdof*(3*i+1)+3*j+0]+=gauss->
weight*Jdet*2*mu_prime*eps1dotdphij*eps2dotdphii;
3116 Ke->
values[numdof*(3*i+1)+3*j+1]+=gauss->
weight*Jdet*2*mu_prime*eps2dotdphij*eps2dotdphii;
3117 Ke->
values[numdof*(3*i+1)+3*j+2]+=gauss->
weight*Jdet*2*mu_prime*eps3dotdphij*eps2dotdphii;
3119 Ke->
values[numdof*(3*i+2)+3*j+0]+=gauss->
weight*Jdet*2*mu_prime*eps1dotdphij*eps3dotdphii;
3120 Ke->
values[numdof*(3*i+2)+3*j+1]+=gauss->
weight*Jdet*2*mu_prime*eps2dotdphij*eps3dotdphii;
3121 Ke->
values[numdof*(3*i+2)+3*j+2]+=gauss->
weight*Jdet*2*mu_prime*eps3dotdphij*eps3dotdphii;
3131 xDelete<IssmDouble>(xyz_list);
3132 xDelete<IssmDouble>(dbasis);
3133 xDelete<int>(cs_list);
◆ CreateKMatrixFS()
◆ CreateKMatrixFSFriction()
Definition at line 3922 of file StressbalanceAnalysis.cpp.
3932 bool mainlyfloating;
3934 int friction_style,point1;
3938 Gauss* gauss = NULL;
3947 int numdof = vnumnodes*dim + pnumnodes;
3951 IssmDouble* vbasis=xNew<IssmDouble>(vnumnodes);
3956 Input2* gllevelset_input = NULL;
3966 element->
GetGroundedPart(&point1,&fraction1,&fraction2,&mainlyfloating);
3975 for(
int ig=gauss->
begin();ig<gauss->end();ig++){
3982 if(gllevelset<0.) alpha2=0.;
3985 if (phi<0.99999999) alpha2=0.;
3987 else _error_(
"friction interpolation "<<
EnumToStringx(friction_style)<<
" not implemented yet");
3994 for(
int i=0;i<vnumnodes;i++){
3995 for(
int j=0;j<vnumnodes;j++){
3996 Ke->
values[(3*i+0)*numdof+3*j+0] += alpha2*gauss->
weight*Jdet*vbasis[i]*vbasis[j];
3997 Ke->
values[(3*i+1)*numdof+3*j+1] += alpha2*gauss->
weight*Jdet*vbasis[i]*vbasis[j];
3998 Ke->
values[(3*i+2)*numdof+3*j+2] += alpha2*gauss->
weight*Jdet*vbasis[i]*vbasis[j];
4004 for(
int i=0;i<vnumnodes;i++){
4005 for(
int j=0;j<vnumnodes;j++){
4006 Ke->
values[(2*i+0)*numdof+2*j+0] += alpha2*gauss->
weight*Jdet*vbasis[i]*vbasis[j];
4007 Ke->
values[(2*i+1)*numdof+2*j+1] += alpha2*gauss->
weight*Jdet*vbasis[i]*vbasis[j];
4018 xDelete<IssmDouble>(xyz_list_base);
4019 xDelete<IssmDouble>(vbasis);
◆ CreateKMatrixFSShelf()
Definition at line 3164 of file StressbalanceAnalysis.cpp.
3169 int approximation,shelf_dampening;
3173 if(shelf_dampening==0)
return NULL;
3176 bool mainlyfloating;
3182 Gauss* gauss = NULL;
3190 int numdof = vnumnodes*dim + pnumnodes;
3194 IssmDouble* vbasis = xNew<IssmDouble>(vnumnodes);
3207 for(
int ig=gauss->
begin();ig<gauss->end();ig++){
3213 if(dim==2) slope2=slope[0]*slope[0];
3214 else if(dim==3) slope2=slope[0]*slope[0]+slope[1]*slope[1];
3215 scalar = rho_water*gravity*sqrt(1+slope2)*gauss->
weight*Jdet*dt;
3216 for(i=0;i<vnumnodes;i++){
3217 for(j=0;j<vnumnodes;j++){
3218 Ke->
values[numdof*((i+1)*dim-1)+(j+1)*dim-1] += scalar*vbasis[i]*vbasis[j];
3227 xDelete<IssmDouble>(xyz_list_base);
3228 xDelete<IssmDouble>(xyz_list);
3229 xDelete<IssmDouble>(vbasis);
◆ CreateKMatrixFSViscous()
Definition at line 3232 of file StressbalanceAnalysis.cpp.
3245 int numdof = vnumnodes*dim + pnumnodes;
3248 int* cs_list = xNew<int>(vnumnodes+pnumnodes);
3249 if(dim==2)
for(i=0;i<vnumnodes;i++) cs_list[i] =
XYEnum;
3250 else for(i=0;i<vnumnodes;i++) cs_list[i] =
XYZEnum;
3251 for(i=0;i<pnumnodes;i++) cs_list[vnumnodes+i] =
PressureEnum;
3255 IssmDouble* vdbasis = xNew<IssmDouble>(dim*vnumnodes);
3256 IssmDouble* pbasis = xNew<IssmDouble>(pnumnodes);
3268 for(
int ig=gauss->
begin();ig<gauss->end();ig++){
3278 for(
int i=0;i<vnumnodes;i++){
3279 for(
int j=0;j<vnumnodes;j++){
3280 Ke->
values[(3*i+0)*numdof+3*j+0] += gauss->
weight*Jdet*viscosity*(2.*vdbasis[0*vnumnodes+j]*vdbasis[0*vnumnodes+i] + vdbasis[1*vnumnodes+j]*vdbasis[1*vnumnodes+i] + vdbasis[2*vnumnodes+j]*vdbasis[2*vnumnodes+i]);
3281 Ke->
values[(3*i+0)*numdof+3*j+1] += gauss->
weight*Jdet*viscosity*(vdbasis[0*vnumnodes+j]*vdbasis[1*vnumnodes+i]);
3282 Ke->
values[(3*i+0)*numdof+3*j+2] += gauss->
weight*Jdet*viscosity*(vdbasis[0*vnumnodes+j]*vdbasis[2*vnumnodes+i]);
3283 Ke->
values[(3*i+1)*numdof+3*j+0] += gauss->
weight*Jdet*viscosity*(vdbasis[1*vnumnodes+j]*vdbasis[0*vnumnodes+i]);
3284 Ke->
values[(3*i+1)*numdof+3*j+1] += gauss->
weight*Jdet*viscosity*(vdbasis[0*vnumnodes+j]*vdbasis[0*vnumnodes+i] + 2.*vdbasis[1*vnumnodes+j]*vdbasis[1*vnumnodes+i] + vdbasis[2*vnumnodes+j]*vdbasis[2*vnumnodes+i]);
3285 Ke->
values[(3*i+1)*numdof+3*j+2] += gauss->
weight*Jdet*viscosity*(vdbasis[1*vnumnodes+j]*vdbasis[2*vnumnodes+i]);
3286 Ke->
values[(3*i+2)*numdof+3*j+0] += gauss->
weight*Jdet*viscosity*(vdbasis[2*vnumnodes+j]*vdbasis[0*vnumnodes+i]);
3287 Ke->
values[(3*i+2)*numdof+3*j+1] += gauss->
weight*Jdet*viscosity*(vdbasis[2*vnumnodes+j]*vdbasis[1*vnumnodes+i]);
3288 Ke->
values[(3*i+2)*numdof+3*j+2] += gauss->
weight*Jdet*viscosity*(vdbasis[0*vnumnodes+j]*vdbasis[0*vnumnodes+i] + vdbasis[1*vnumnodes+j]*vdbasis[1*vnumnodes+i] + 2.*vdbasis[2*vnumnodes+j]*vdbasis[2*vnumnodes+i]);
3290 for(
int k=0;k<pnumnodes;k++){
3291 Ke->
values[(3*i+0)*numdof+3*vnumnodes+k] += gauss->
weight*Jdet*FSreconditioning*(-pbasis[k]*vdbasis[0*vnumnodes+i]);
3292 Ke->
values[(3*i+1)*numdof+3*vnumnodes+k] += gauss->
weight*Jdet*FSreconditioning*(-pbasis[k]*vdbasis[1*vnumnodes+i]);
3293 Ke->
values[(3*i+2)*numdof+3*vnumnodes+k] += gauss->
weight*Jdet*FSreconditioning*(-pbasis[k]*vdbasis[2*vnumnodes+i]);
3297 for(
int k=0;k<pnumnodes;k++){
3298 for(
int j=0;j<vnumnodes;j++){
3299 Ke->
values[(3*vnumnodes+k)*numdof+3*j+0] += gauss->
weight*Jdet*(-FSreconditioning*vdbasis[0*vnumnodes+j]*pbasis[k]);
3300 Ke->
values[(3*vnumnodes+k)*numdof+3*j+1] += gauss->
weight*Jdet*(-FSreconditioning*vdbasis[1*vnumnodes+j]*pbasis[k]);
3301 Ke->
values[(3*vnumnodes+k)*numdof+3*j+2] += gauss->
weight*Jdet*(-FSreconditioning*vdbasis[2*vnumnodes+j]*pbasis[k]);
3307 for(
int i=0;i<vnumnodes;i++){
3308 for(
int j=0;j<vnumnodes;j++){
3309 Ke->
values[(2*i+0)*numdof+2*j+0] += gauss->
weight*Jdet*viscosity*(2.*vdbasis[0*vnumnodes+j]*vdbasis[0*vnumnodes+i] + vdbasis[1*vnumnodes+j]*vdbasis[1*vnumnodes+i]);
3310 Ke->
values[(2*i+0)*numdof+2*j+1] += gauss->
weight*Jdet*viscosity*(vdbasis[0*vnumnodes+j]*vdbasis[1*vnumnodes+i]);
3311 Ke->
values[(2*i+1)*numdof+2*j+0] += gauss->
weight*Jdet*viscosity*(vdbasis[1*vnumnodes+j]*vdbasis[0*vnumnodes+i]);
3312 Ke->
values[(2*i+1)*numdof+2*j+1] += gauss->
weight*Jdet*viscosity*(vdbasis[0*vnumnodes+j]*vdbasis[0*vnumnodes+i] + 2.*vdbasis[1*vnumnodes+j]*vdbasis[1*vnumnodes+i]);
3314 for(
int k=0;k<pnumnodes;k++){
3315 Ke->
values[(2*i+0)*numdof+2*vnumnodes+k] += gauss->
weight*Jdet*FSreconditioning*(-pbasis[k]*vdbasis[0*vnumnodes+i]);
3316 Ke->
values[(2*i+1)*numdof+2*vnumnodes+k] += gauss->
weight*Jdet*FSreconditioning*(-pbasis[k]*vdbasis[1*vnumnodes+i]);
3320 for(
int k=0;k<pnumnodes;k++){
3321 for(
int j=0;j<vnumnodes;j++){
3322 Ke->
values[(2*vnumnodes+k)*numdof+2*j+0] += gauss->
weight*Jdet*(-FSreconditioning*vdbasis[0*vnumnodes+j]*pbasis[k]);
3323 Ke->
values[(2*vnumnodes+k)*numdof+2*j+1] += gauss->
weight*Jdet*(-FSreconditioning*vdbasis[1*vnumnodes+j]*pbasis[k]);
3334 xDelete<IssmDouble>(xyz_list);
3335 xDelete<IssmDouble>(pbasis);
3336 xDelete<IssmDouble>(vdbasis);
3337 xDelete<int>(cs_list);
◆ CreateKMatrixFSViscousLA()
Definition at line 3446 of file StressbalanceAnalysis.cpp.
3458 if(dim==2) epssize = 3;
3465 int numdof = vnumnodes*dim;
3466 int pnumdof = pnumnodes;
3467 int lnumdof = lnumnodes;
3470 int* cs_list = xNew<int>(vnumnodes);
3471 if(dim==2)
for(i=0;i<vnumnodes;i++) cs_list[i] =
XYEnum;
3472 else for(i=0;i<vnumnodes;i++) cs_list[i] =
XYZEnum;
3476 IssmDouble* B = xNew<IssmDouble>(epssize*numdof);
3477 IssmDouble* Bprime = xNew<IssmDouble>(epssize*numdof);
3478 IssmDouble* BtBUzawa = xNewZeroInit<IssmDouble>(numdof*pnumdof);
3480 IssmDouble* BprimeU = xNew<IssmDouble>(numdof);
3481 IssmDouble* D = xNewZeroInit<IssmDouble>(epssize*epssize);
3492 for(
int ig=gauss->
begin();ig<gauss->end();ig++){
3496 this->
GetBFSvel(B,element,dim,xyz_list,gauss);
3500 for(i=0;i<epssize;i++) D[i*epssize+i] = 2*viscosity*gauss->
weight*Jdet;
3503 D,epssize,epssize,0,
3504 Bprime,epssize,numdof,0,
3510 DU = gauss->
weight*Jdet*sqrt(r);
3520 BtBUzawa,pnumdof,numdof,0,
3526 element->
NormalBase(&normal[0],xyz_list_base);
3528 IssmDouble* Dlambda = xNewZeroInit<IssmDouble>(dim*dim);
3529 IssmDouble* C = xNewZeroInit<IssmDouble>(dim*lnumdof);
3530 IssmDouble* Cprime = xNewZeroInit<IssmDouble>(dim*numdof);
3531 IssmDouble* CtCUzawa = xNewZeroInit<IssmDouble>(numdof*lnumdof);
3535 for(
int ig=gauss->
begin();ig<gauss->end();ig++){
3539 this->
GetCFS(C,element,dim,xyz_list,gauss);
3540 this->
GetCFSprime(Cprime,element,dim,xyz_list,gauss);
3541 for(i=0;i<dim;i++) Dlambda[i*dim+i] = gauss->
weight*Jdet*sqrt(normal[i]*normal[i])*sqrt(rl);
3544 Cprime,dim,numdof,0,
3550 CtCUzawa,lnumdof,numdof,0,
3554 xDelete<IssmDouble>(Dlambda);
3555 xDelete<IssmDouble>(C);
3556 xDelete<IssmDouble>(Cprime);
3557 xDelete<IssmDouble>(CtCUzawa);
3558 xDelete<IssmDouble>(xyz_list_base);
3566 xDelete<IssmDouble>(xyz_list);
3567 xDelete<IssmDouble>(D);
3568 xDelete<IssmDouble>(Bprime);
3569 xDelete<IssmDouble>(B);
3570 xDelete<IssmDouble>(BprimeU);
3571 xDelete<IssmDouble>(BU);
3572 xDelete<IssmDouble>(BtBUzawa);
3573 xDelete<int>(cs_list);
◆ CreateKMatrixFSViscousXTH()
Definition at line 3576 of file StressbalanceAnalysis.cpp.
3586 if(dim==2) epssize = 3;
3592 int numdof = vnumnodes*dim + pnumnodes;
3593 int bsize = epssize + 2;
3596 int* cs_list = xNew<int>(vnumnodes+pnumnodes);
3597 if(dim==2)
for(i=0;i<vnumnodes;i++) cs_list[i] =
XYEnum;
3598 else for(i=0;i<vnumnodes;i++) cs_list[i] =
XYZEnum;
3599 for(i=0;i<pnumnodes;i++) cs_list[vnumnodes+i] =
PressureEnum;
3603 IssmDouble* B = xNew<IssmDouble>(bsize*numdof);
3604 IssmDouble* Bprime = xNew<IssmDouble>(bsize*numdof);
3605 IssmDouble* D = xNewZeroInit<IssmDouble>(bsize*bsize);
3617 for(
int ig=gauss->
begin();ig<gauss->end();ig++){
3621 this->
GetBFS(B,element,dim,xyz_list,gauss);
3622 this->
GetBFSprime(Bprime,element,dim,xyz_list,gauss);
3624 for(i=0;i<epssize;i++) D[i*bsize+i] = + r*gauss->
weight*Jdet;
3625 for(i=epssize;i<bsize;i++) D[i*bsize+i] = - FSreconditioning*gauss->
weight*Jdet;
3629 Bprime,bsize,numdof,0,
3638 xDelete<IssmDouble>(xyz_list);
3639 xDelete<IssmDouble>(D);
3640 xDelete<IssmDouble>(Bprime);
3641 xDelete<IssmDouble>(B);
3642 xDelete<int>(cs_list);
◆ CreatePressureMassMatrix()
Definition at line 3340 of file StressbalanceAnalysis.cpp.
3353 int numdof = vnumnodes*dim + pnumnodes;
3357 IssmDouble* pbasis = xNew<IssmDouble>(pnumnodes);
3365 for(
int ig=gauss->
begin();ig<gauss->end();ig++){
3371 if(dim==3 || dim==2){
3373 for(
int k=0;k<pnumnodes;k++){
3374 for(
int j=0;j<pnumnodes;j++){
3375 Ke->
values[(dim*vnumnodes+k)*numdof+dim*vnumnodes+j] += gauss->
weight*Jdet*(pbasis[j]*pbasis[k]);
3385 xDelete<IssmDouble>(xyz_list);
3386 xDelete<IssmDouble>(pbasis);
◆ CreateSchurPrecondMatrix()
Definition at line 3389 of file StressbalanceAnalysis.cpp.
3402 int numdof = vnumnodes*dim + pnumnodes;
3406 IssmDouble* pbasis = xNew<IssmDouble>(pnumnodes);
3419 for(
int ig=gauss->
begin();ig<gauss->end();ig++){
3427 if(dim==3 || dim==2){
3429 for(
int k=0;k<pnumnodes;k++){
3430 for(
int j=0;j<pnumnodes;j++){
3431 Ke->
values[(dim*vnumnodes+k)*numdof+dim*vnumnodes+j] += gauss->
weight*1./viscosity*Jdet*(pbasis[j]*pbasis[k]);
3441 xDelete<IssmDouble>(xyz_list);
3442 xDelete<IssmDouble>(pbasis);
◆ CreatePVectorFS()
◆ CreatePVectorFSFriction()
◆ CreatePVectorFSFront()
Definition at line 4140 of file StressbalanceAnalysis.cpp.
4151 Gauss *gauss = NULL;
4162 int* cs_list = xNew<int>(vnumnodes+pnumnodes);
4163 if(dim==2)
for(i=0;i<vnumnodes;i++) cs_list[i] =
XYEnum;
4164 else for(i=0;i<vnumnodes;i++) cs_list[i] =
XYZEnum;
4165 for(i=0;i<pnumnodes;i++) cs_list[vnumnodes+i] =
PressureEnum;
4169 IssmDouble* vbasis = xNew<IssmDouble>(vnumnodes);
4181 IssmDouble zmax=xyz_list[0*3+(dim-1)];
for(
int i=1;i<numvertices;i++)
if(xyz_list[i*3+(dim-1)]>zmax) zmax=xyz_list[i*3+(dim-1)];
4182 IssmDouble zmin=xyz_list[0*3+(dim-1)];
for(
int i=1;i<numvertices;i++)
if(xyz_list[i*3+(dim-1)]<zmin) zmin=xyz_list[i*3+(dim-1)];
4183 if(zmax>0. && zmin<0.) gauss=element->
NewGauss(xyz_list,xyz_list_front,3,30);
4184 else gauss=element->
NewGauss(xyz_list,xyz_list_front,3,3);
4187 for(
int ig=gauss->
begin();ig<gauss->
end();ig++){
4194 if(dim==3) z=element->
GetZcoord(xyz_list,gauss);
4195 else z=element->
GetYcoord(xyz_list,gauss);
4196 pressure = rho_water*gravity*
min(0.,z-sealevel);
4198 for (
int i=0;i<vnumnodes;i++){
4199 pe->
values[dim*i+0]+= pressure*Jdet*gauss->
weight*normal[0]*vbasis[i];
4200 pe->
values[dim*i+1]+= pressure*Jdet*gauss->
weight*normal[1]*vbasis[i];
4201 if(dim==3) pe->
values[dim*i+2]+= pressure*Jdet*gauss->
weight*normal[2]*vbasis[i];
4210 xDelete<int>(cs_list);
4211 xDelete<IssmDouble>(vbasis);
4212 xDelete<IssmDouble>(xyz_list);
4213 xDelete<IssmDouble>(xyz_list_front);
◆ CreatePVectorFSShelf()
Definition at line 4216 of file StressbalanceAnalysis.cpp.
4233 int* cs_list = xNew<int>(vnumnodes+pnumnodes);
4234 if(dim==2)
for(i=0;i<vnumnodes;i++) cs_list[i] =
XYEnum;
4235 else for(i=0;i<vnumnodes;i++) cs_list[i] =
XYZEnum;
4236 for(i=0;i<pnumnodes;i++) cs_list[vnumnodes+i] =
PressureEnum;
4240 IssmDouble* vbasis = xNew<IssmDouble>(vnumnodes);
4250 for(
int ig=gauss->
begin();ig<gauss->end();ig++){
4257 water_pressure=gravity*rho_water*base;
4259 for(i=0;i<vnumnodes;i++){
4260 pe->
values[i*dim+(dim-1)]+=-water_pressure*gauss->
weight*Jdet*vbasis[i];
4265 int shelf_dampening;
4267 if(shelf_dampening) {
4271 for(
int ig=gauss->
begin();ig<gauss->end();ig++){
4276 for(i=0;i<vnumnodes;i++){
4277 pe->
values[i*dim+(dim-1)] += -dt*rho_water*gravity*mb*gauss->
weight*Jdet*vbasis[i];
4286 xDelete<int>(cs_list);
4287 xDelete<IssmDouble>(vbasis);
4288 xDelete<IssmDouble>(xyz_list_base);
◆ CreatePVectorFSStress()
◆ CreatePVectorFSViscous()
Definition at line 4070 of file StressbalanceAnalysis.cpp.
4084 int* cs_list = xNew<int>(vnumnodes+pnumnodes);
4085 if(dim==2)
for(i=0;i<vnumnodes;i++) cs_list[i] =
XYEnum;
4086 else for(i=0;i<vnumnodes;i++) cs_list[i] =
XYZEnum;
4087 for(i=0;i<pnumnodes;i++) cs_list[vnumnodes+i] =
PressureEnum;
4091 IssmDouble* vbasis = xNew<IssmDouble>(vnumnodes);
4099 Input2* loadingforcez_input=NULL;
4106 for(
int ig=gauss->
begin();ig<gauss->end();ig++){
4114 if(dim==3) loadingforcez_input->
GetInputValue(&forcez,gauss);
4116 for(i=0;i<vnumnodes;i++){
4117 pe->
values[i*dim+0] += +rho_ice*forcex *Jdet*gauss->
weight*vbasis[i];
4118 pe->
values[i*dim+1] += +rho_ice*forcey *Jdet*gauss->
weight*vbasis[i];
4120 pe->
values[i*dim+2] += +rho_ice*forcez*Jdet*gauss->
weight*vbasis[i];
4121 pe->
values[i*dim+2] += -rho_ice*gravity*Jdet*gauss->
weight*vbasis[i];
4124 pe->
values[i*dim+1] += -rho_ice*gravity*Jdet*gauss->
weight*vbasis[i];
4134 xDelete<int>(cs_list);
4135 xDelete<IssmDouble>(vbasis);
4136 xDelete<IssmDouble>(xyz_list);
◆ CreatePVectorFSViscousLA()
Definition at line 4291 of file StressbalanceAnalysis.cpp.
4298 Gauss* gauss = NULL;
4307 int* cs_list = xNew<int>(numnodes);
4308 if(dim==2)
for(i=0;i<numnodes;i++) cs_list[i] =
XYEnum;
4309 else for(i=0;i<numnodes;i++) cs_list[i] =
XYZEnum;
4313 IssmDouble* dbasis = xNew<IssmDouble>(3*numnodes);
4323 for(
int ig=gauss->
begin();ig<gauss->end();ig++){
4330 for(i=0;i<numnodes;i++){
4331 pe->
values[i*dim+0] += pressure*gauss->
weight*Jdet*dbasis[0*numnodes+i];
4332 pe->
values[i*dim+1] += pressure*gauss->
weight*Jdet*dbasis[1*numnodes+i];
4333 if(dim==3) pe->
values[i*dim+2]+= pressure*gauss->
weight*Jdet*dbasis[2*numnodes+i];
4339 IssmDouble* vbasis = xNew<IssmDouble>(numnodes);
4342 element->
NormalBase(&bed_normal[0],xyz_list_base);
4346 for(
int ig=gauss->
begin();ig<gauss->end();ig++){
4353 for(i=0;i<numnodes;i++){
4354 pe->
values[i*dim+0] += + sigmann*bed_normal[0]*gauss->
weight*Jdet*vbasis[i];
4355 pe->
values[i*dim+1] += + sigmann*bed_normal[1]*gauss->
weight*Jdet*vbasis[i];
4356 if(dim==3) pe->
values[i*dim+2] += + sigmann*bed_normal[2]*gauss->
weight*Jdet*vbasis[i];
4359 xDelete<IssmDouble>(xyz_list_base);
4360 xDelete<IssmDouble>(vbasis);
4368 xDelete<int>(cs_list);
4369 xDelete<IssmDouble>(xyz_list);
4370 xDelete<IssmDouble>(dbasis);
◆ CreatePVectorFSViscousXTH()
Definition at line 4373 of file StressbalanceAnalysis.cpp.
4377 IssmDouble epsxx,epsyy,epszz,epsxy,epsxz,epsyz;
4378 IssmDouble sigmapxx,sigmapyy,sigmapzz,sigmapxy,sigmapxz,sigmapyz;
4380 Gauss* gauss = NULL;
4384 if(dim==2) tausize = 3;
4393 int* cs_list = xNew<int>(vnumnodes+pnumnodes);
4394 if(dim==2)
for(i=0;i<vnumnodes;i++) cs_list[i] =
XYEnum;
4395 else for(i=0;i<vnumnodes;i++) cs_list[i] =
XYZEnum;
4396 for(i=0;i<pnumnodes;i++) cs_list[vnumnodes+i] =
PressureEnum;
4400 IssmDouble* Dstar = xNewZeroInit<IssmDouble>((dim*vnumnodes)*(tausize*tnumnodes));
4401 IssmDouble* tau = xNew<IssmDouble>(tausize*tnumnodes);
4402 IssmDouble* d = xNew<IssmDouble>(tausize*tnumnodes);
4403 IssmDouble* vdbasis = xNew<IssmDouble>(dim*vnumnodes);
4404 IssmDouble* tbasis = xNew<IssmDouble>(tnumnodes);
4405 IssmDouble* D = xNewZeroInit<IssmDouble>(tausize*tnumnodes*tausize*tnumnodes);
4419 Input2* sigmapzz_input=NULL;
Input2* sigmapxz_input=NULL;
Input2* sigmapyz_input=NULL;
4430 for(
int i=0;i<tnumnodes;i++){
4437 d[i*tausize+0]=epsxx; tau[i*tausize+0]=sigmapxx;
4438 d[i*tausize+1]=epsyy; tau[i*tausize+1]=sigmapyy;
4439 d[i*tausize+2]=epsxy; tau[i*tausize+2]=sigmapxy;
4445 d[i*tausize+0]=epsxx; tau[i*tausize+0]=sigmapxx;
4446 d[i*tausize+1]=epsyy; tau[i*tausize+1]=sigmapyy;
4447 d[i*tausize+2]=epszz; tau[i*tausize+2]=sigmapzz;
4448 d[i*tausize+3]=epsxy; tau[i*tausize+3]=sigmapxy;
4449 d[i*tausize+4]=epsxz; tau[i*tausize+4]=sigmapxz;
4450 d[i*tausize+5]=epsyz; tau[i*tausize+5]=sigmapyz;
4457 for(
int ig=gauss->
begin();ig<gauss->end();ig++){
4487 for(
int i=0;i<vnumnodes;i++){
4488 for(
int j=0;j<tnumnodes;j++){
4489 Dstar[(i*dim+0)*tausize*tnumnodes + j*tausize+0] += gauss->
weight*Jdet*tbasis[j]*vdbasis[0*vnumnodes+i];
4490 Dstar[(i*dim+0)*tausize*tnumnodes + j*tausize+2] += gauss->
weight*Jdet*tbasis[j]*vdbasis[1*vnumnodes+i];
4492 Dstar[(i*dim+1)*tausize*tnumnodes + j*tausize+1] += gauss->
weight*Jdet*tbasis[j]*vdbasis[1*vnumnodes+i];
4493 Dstar[(i*dim+1)*tausize*tnumnodes + j*tausize+2] += gauss->
weight*Jdet*tbasis[j]*vdbasis[0*vnumnodes+i];
4498 for(
int i=0;i<vnumnodes;i++){
4499 for(
int j=0;j<tnumnodes;j++){
4500 Dstar[(i*dim+0)*tausize*tnumnodes + j*tausize+0] += gauss->
weight*Jdet*tbasis[j]*vdbasis[0*vnumnodes+i];
4501 Dstar[(i*dim+0)*tausize*tnumnodes + j*tausize+3] += gauss->
weight*Jdet*tbasis[j]*vdbasis[1*vnumnodes+i];
4502 Dstar[(i*dim+0)*tausize*tnumnodes + j*tausize+4] += gauss->
weight*Jdet*tbasis[j]*vdbasis[2*vnumnodes+i];
4504 Dstar[(i*dim+1)*tausize*tnumnodes + j*tausize+1] += gauss->
weight*Jdet*tbasis[j]*vdbasis[1*vnumnodes+i];
4505 Dstar[(i*dim+1)*tausize*tnumnodes + j*tausize+3] += gauss->
weight*Jdet*tbasis[j]*vdbasis[0*vnumnodes+i];
4506 Dstar[(i*dim+1)*tausize*tnumnodes + j*tausize+5] += gauss->
weight*Jdet*tbasis[j]*vdbasis[2*vnumnodes+i];
4508 Dstar[(i*dim+2)*tausize*tnumnodes + j*tausize+2] += gauss->
weight*Jdet*tbasis[j]*vdbasis[2*vnumnodes+i];
4509 Dstar[(i*dim+2)*tausize*tnumnodes + j*tausize+4] += gauss->
weight*Jdet*tbasis[j]*vdbasis[0*vnumnodes+i];
4510 Dstar[(i*dim+2)*tausize*tnumnodes + j*tausize+5] += gauss->
weight*Jdet*tbasis[j]*vdbasis[1*vnumnodes+i];
4517 for(i=0;i<tausize*tnumnodes;i++) D[i*(tausize*tnumnodes)+i] = -1.;
4519 D,tausize*tnumnodes,tausize*tnumnodes,0,
4520 tau,tausize*tnumnodes,1,0,
4524 for(i=0;i<tausize*tnumnodes;i++) D[i*(tausize*tnumnodes)+i] = +r;
4526 D,tausize*tnumnodes,tausize*tnumnodes,0,
4527 d,tausize*tnumnodes,1,0,
4535 xDelete<int>(cs_list);
4536 xDelete<IssmDouble>(xyz_list);
4537 xDelete<IssmDouble>(Dstar);
4538 xDelete<IssmDouble>(d);
4539 xDelete<IssmDouble>(D);
4540 xDelete<IssmDouble>(tau);
4541 xDelete<IssmDouble>(vdbasis);
4542 xDelete<IssmDouble>(tbasis);
◆ GetBFS()
Definition at line 4545 of file StressbalanceAnalysis.cpp.
4587 IssmDouble* vdbasis=xNew<IssmDouble>(dim*vnumnodes);
4588 IssmDouble* pbasis =xNew<IssmDouble>(pnumnodes);
4594 for(
int i=0;i<vnumnodes;i++){
4595 B[(dim*vnumnodes+pnumnodes)*0+dim*i+0] = vdbasis[0*vnumnodes+i];
4596 B[(dim*vnumnodes+pnumnodes)*0+dim*i+1] = 0.;
4597 B[(dim*vnumnodes+pnumnodes)*1+dim*i+0] = 0.;
4598 B[(dim*vnumnodes+pnumnodes)*1+dim*i+1] = vdbasis[1*vnumnodes+i];
4599 B[(dim*vnumnodes+pnumnodes)*2+dim*i+0] = .5*vdbasis[1*vnumnodes+i];
4600 B[(dim*vnumnodes+pnumnodes)*2+dim*i+1] = .5*vdbasis[0*vnumnodes+i];
4601 B[(dim*vnumnodes+pnumnodes)*3+dim*i+0] = 0.;
4602 B[(dim*vnumnodes+pnumnodes)*3+dim*i+1] = 0.;
4603 B[(dim*vnumnodes+pnumnodes)*4+dim*i+0] = vdbasis[0*vnumnodes+i];
4604 B[(dim*vnumnodes+pnumnodes)*4+dim*i+1] = vdbasis[1*vnumnodes+i];
4606 for(
int i=0;i<pnumnodes;i++){
4607 B[(dim*vnumnodes+pnumnodes)*0+(dim*vnumnodes)+i] = 0.;
4608 B[(dim*vnumnodes+pnumnodes)*1+(dim*vnumnodes)+i] = 0.;
4609 B[(dim*vnumnodes+pnumnodes)*2+(dim*vnumnodes)+i] = 0.;
4610 B[(dim*vnumnodes+pnumnodes)*3+(dim*vnumnodes)+i] = pbasis[i];
4611 B[(dim*vnumnodes+pnumnodes)*4+(dim*vnumnodes)+i] = 0.;
4615 for(
int i=0;i<vnumnodes;i++){
4616 B[(dim*vnumnodes+pnumnodes)*0+dim*i+0] = vdbasis[0*vnumnodes+i];
4617 B[(dim*vnumnodes+pnumnodes)*0+dim*i+1] = 0.;
4618 B[(dim*vnumnodes+pnumnodes)*0+dim*i+2] = 0.;
4619 B[(dim*vnumnodes+pnumnodes)*1+dim*i+0] = 0.;
4620 B[(dim*vnumnodes+pnumnodes)*1+dim*i+1] = vdbasis[1*vnumnodes+i];
4621 B[(dim*vnumnodes+pnumnodes)*1+dim*i+2] = 0.;
4622 B[(dim*vnumnodes+pnumnodes)*2+dim*i+0] = 0.;
4623 B[(dim*vnumnodes+pnumnodes)*2+dim*i+1] = 0.;
4624 B[(dim*vnumnodes+pnumnodes)*2+dim*i+2] = vdbasis[2*vnumnodes+i];
4625 B[(dim*vnumnodes+pnumnodes)*3+dim*i+0] = .5*vdbasis[1*vnumnodes+i];
4626 B[(dim*vnumnodes+pnumnodes)*3+dim*i+1] = .5*vdbasis[0*vnumnodes+i];
4627 B[(dim*vnumnodes+pnumnodes)*3+dim*i+2] = 0.;
4628 B[(dim*vnumnodes+pnumnodes)*4+dim*i+0] = .5*vdbasis[2*vnumnodes+i];
4629 B[(dim*vnumnodes+pnumnodes)*4+dim*i+1] = 0.;
4630 B[(dim*vnumnodes+pnumnodes)*4+dim*i+2] = .5*vdbasis[0*vnumnodes+i];
4631 B[(dim*vnumnodes+pnumnodes)*5+dim*i+0] = 0.;
4632 B[(dim*vnumnodes+pnumnodes)*5+dim*i+1] = .5*vdbasis[2*vnumnodes+i];
4633 B[(dim*vnumnodes+pnumnodes)*5+dim*i+2] = .5*vdbasis[1*vnumnodes+i];
4634 B[(dim*vnumnodes+pnumnodes)*6+dim*i+0] = 0.;
4635 B[(dim*vnumnodes+pnumnodes)*6+dim*i+1] = 0.;
4636 B[(dim*vnumnodes+pnumnodes)*6+dim*i+2] = 0.;
4637 B[(dim*vnumnodes+pnumnodes)*7+dim*i+0] = vdbasis[0*vnumnodes+i];
4638 B[(dim*vnumnodes+pnumnodes)*7+dim*i+1] = vdbasis[1*vnumnodes+i];
4639 B[(dim*vnumnodes+pnumnodes)*7+dim*i+2] = vdbasis[2*vnumnodes+i];
4641 for(
int i=0;i<pnumnodes;i++){
4642 B[(dim*vnumnodes+pnumnodes)*0+(dim*vnumnodes)+i] = 0.;
4643 B[(dim*vnumnodes+pnumnodes)*1+(dim*vnumnodes)+i] = 0.;
4644 B[(dim*vnumnodes+pnumnodes)*2+(dim*vnumnodes)+i] = 0.;
4645 B[(dim*vnumnodes+pnumnodes)*3+(dim*vnumnodes)+i] = 0.;
4646 B[(dim*vnumnodes+pnumnodes)*4+(dim*vnumnodes)+i] = 0.;
4647 B[(dim*vnumnodes+pnumnodes)*5+(dim*vnumnodes)+i] = 0.;
4648 B[(dim*vnumnodes+pnumnodes)*6+(dim*vnumnodes)+i] = pbasis[i];
4649 B[(dim*vnumnodes+pnumnodes)*7+(dim*vnumnodes)+i] = 0.;
4654 xDelete<IssmDouble>(vdbasis);
4655 xDelete<IssmDouble>(pbasis);
◆ GetBFSFriction()
Definition at line 4657 of file StressbalanceAnalysis.cpp.
4671 int pnumdof = pnumnodes;
4672 int vnumdof = vnumnodes*dim;
4675 IssmDouble* vbasis=xNew<IssmDouble>(vnumnodes);
4680 for(
int i=0;i<vnumnodes;i++){
4681 B[(vnumdof+pnumdof)*0+3*i+0] = vbasis[i];
4682 B[(vnumdof+pnumdof)*0+3*i+1] = 0.;
4683 B[(vnumdof+pnumdof)*0+3*i+2] = 0.;
4685 B[(vnumdof+pnumdof)*1+3*i+0] = 0.;
4686 B[(vnumdof+pnumdof)*1+3*i+1] = vbasis[i];
4687 B[(vnumdof+pnumdof)*1+3*i+2] = 0.;
4689 B[(vnumdof+pnumdof)*2+3*i+0] = 0.;
4690 B[(vnumdof+pnumdof)*2+3*i+1] = 0.;
4691 B[(vnumdof+pnumdof)*2+3*i+2] = vbasis[i];
4693 for(
int i=0;i<pnumnodes;i++){
4694 B[(vnumdof+pnumdof)*0+i+vnumdof+0] = 0.;
4695 B[(vnumdof+pnumdof)*1+i+vnumdof+0] = 0.;
4696 B[(vnumdof+pnumdof)*2+i+vnumdof+0] = 0.;
4700 for(
int i=0;i<vnumnodes;i++){
4701 B[(vnumdof+pnumdof)*0+2*i+0] = vbasis[i];
4702 B[(vnumdof+pnumdof)*0+2*i+1] = 0.;
4704 B[(vnumdof+pnumdof)*1+2*i+0] = 0.;
4705 B[(vnumdof+pnumdof)*1+2*i+1] = vbasis[i];
4708 for(
int i=0;i<pnumnodes;i++){
4709 B[(vnumdof+pnumdof)*0+i+vnumdof+0] = 0.;
4710 B[(vnumdof+pnumdof)*1+i+vnumdof+0] = 0.;
4715 xDelete<IssmDouble>(vbasis);
◆ GetBFSprime()
Definition at line 4717 of file StressbalanceAnalysis.cpp.
4760 IssmDouble* vdbasis=xNew<IssmDouble>(dim*vnumnodes);
4761 IssmDouble* pbasis =xNew<IssmDouble>(pnumnodes);
4767 for(
int i=0;i<vnumnodes;i++){
4768 Bprime[(dim*vnumnodes+pnumnodes)*0+dim*i+0] = vdbasis[0*vnumnodes+i];
4769 Bprime[(dim*vnumnodes+pnumnodes)*0+dim*i+1] = 0.;
4770 Bprime[(dim*vnumnodes+pnumnodes)*1+dim*i+0] = 0.;
4771 Bprime[(dim*vnumnodes+pnumnodes)*1+dim*i+1] = vdbasis[1*vnumnodes+i];
4772 Bprime[(dim*vnumnodes+pnumnodes)*2+dim*i+0] = vdbasis[1*vnumnodes+i];
4773 Bprime[(dim*vnumnodes+pnumnodes)*2+dim*i+1] = vdbasis[0*vnumnodes+i];
4774 Bprime[(dim*vnumnodes+pnumnodes)*3+dim*i+0] = vdbasis[0*vnumnodes+i];
4775 Bprime[(dim*vnumnodes+pnumnodes)*3+dim*i+1] = vdbasis[1*vnumnodes+i];
4776 Bprime[(dim*vnumnodes+pnumnodes)*4+dim*i+0] = 0.;
4777 Bprime[(dim*vnumnodes+pnumnodes)*4+dim*i+1] = 0.;
4779 for(
int i=0;i<pnumnodes;i++){
4780 Bprime[(dim*vnumnodes+pnumnodes)*0+(dim*vnumnodes)+i] = 0.;
4781 Bprime[(dim*vnumnodes+pnumnodes)*1+(dim*vnumnodes)+i] = 0.;
4782 Bprime[(dim*vnumnodes+pnumnodes)*2+(dim*vnumnodes)+i] = 0.;
4783 Bprime[(dim*vnumnodes+pnumnodes)*3+(dim*vnumnodes)+i] = 0.;
4784 Bprime[(dim*vnumnodes+pnumnodes)*4+(dim*vnumnodes)+i] = pbasis[i];
4788 for(
int i=0;i<vnumnodes;i++){
4789 Bprime[(dim*vnumnodes+pnumnodes)*0+dim*i+0] = vdbasis[0*vnumnodes+i];
4790 Bprime[(dim*vnumnodes+pnumnodes)*0+dim*i+1] = 0.;
4791 Bprime[(dim*vnumnodes+pnumnodes)*0+dim*i+2] = 0.;
4792 Bprime[(dim*vnumnodes+pnumnodes)*1+dim*i+0] = 0.;
4793 Bprime[(dim*vnumnodes+pnumnodes)*1+dim*i+1] = vdbasis[1*vnumnodes+i];
4794 Bprime[(dim*vnumnodes+pnumnodes)*1+dim*i+2] = 0.;
4795 Bprime[(dim*vnumnodes+pnumnodes)*2+dim*i+0] = 0.;
4796 Bprime[(dim*vnumnodes+pnumnodes)*2+dim*i+1] = 0.;
4797 Bprime[(dim*vnumnodes+pnumnodes)*2+dim*i+2] = vdbasis[2*vnumnodes+i];
4798 Bprime[(dim*vnumnodes+pnumnodes)*3+dim*i+0] = vdbasis[1*vnumnodes+i];
4799 Bprime[(dim*vnumnodes+pnumnodes)*3+dim*i+1] = vdbasis[0*vnumnodes+i];
4800 Bprime[(dim*vnumnodes+pnumnodes)*3+dim*i+2] = 0.;
4801 Bprime[(dim*vnumnodes+pnumnodes)*4+dim*i+0] = vdbasis[2*vnumnodes+i];
4802 Bprime[(dim*vnumnodes+pnumnodes)*4+dim*i+1] = 0.;
4803 Bprime[(dim*vnumnodes+pnumnodes)*4+dim*i+2] = vdbasis[0*vnumnodes+i];
4804 Bprime[(dim*vnumnodes+pnumnodes)*5+dim*i+0] = 0.;
4805 Bprime[(dim*vnumnodes+pnumnodes)*5+dim*i+1] = vdbasis[2*vnumnodes+i];
4806 Bprime[(dim*vnumnodes+pnumnodes)*5+dim*i+2] = vdbasis[1*vnumnodes+i];
4807 Bprime[(dim*vnumnodes+pnumnodes)*6+dim*i+0] = vdbasis[0*vnumnodes+i];
4808 Bprime[(dim*vnumnodes+pnumnodes)*6+dim*i+1] = vdbasis[1*vnumnodes+i];
4809 Bprime[(dim*vnumnodes+pnumnodes)*6+dim*i+2] = vdbasis[2*vnumnodes+i];
4810 Bprime[(dim*vnumnodes+pnumnodes)*7+dim*i+0] = 0.;
4811 Bprime[(dim*vnumnodes+pnumnodes)*7+dim*i+1] = 0.;
4812 Bprime[(dim*vnumnodes+pnumnodes)*7+dim*i+2] = 0.;
4814 for(
int i=0;i<pnumnodes;i++){
4815 Bprime[(dim*vnumnodes+pnumnodes)*0+(dim*vnumnodes)+i] = 0.;
4816 Bprime[(dim*vnumnodes+pnumnodes)*1+(dim*vnumnodes)+i] = 0.;
4817 Bprime[(dim*vnumnodes+pnumnodes)*2+(dim*vnumnodes)+i] = 0.;
4818 Bprime[(dim*vnumnodes+pnumnodes)*3+(dim*vnumnodes)+i] = 0.;
4819 Bprime[(dim*vnumnodes+pnumnodes)*4+(dim*vnumnodes)+i] = 0.;
4820 Bprime[(dim*vnumnodes+pnumnodes)*5+(dim*vnumnodes)+i] = 0.;
4821 Bprime[(dim*vnumnodes+pnumnodes)*6+(dim*vnumnodes)+i] = 0.;
4822 Bprime[(dim*vnumnodes+pnumnodes)*7+(dim*vnumnodes)+i] = pbasis[i];
4827 xDelete<IssmDouble>(vdbasis);
4828 xDelete<IssmDouble>(pbasis);
◆ GetBFSprimeUzawa()
Definition at line 4830 of file StressbalanceAnalysis.cpp.
4845 IssmDouble* vdbasis=xNew<IssmDouble>(dim*vnumnodes);
4850 for(
int i=0;i<vnumnodes;i++){
4851 Bprime[dim*i+0] = vdbasis[0*vnumnodes+i];
4852 Bprime[dim*i+1] = vdbasis[1*vnumnodes+i];
4856 for(
int i=0;i<vnumnodes;i++){
4857 Bprime[dim*i+0] = vdbasis[0*vnumnodes+i];
4858 Bprime[dim*i+1] = vdbasis[1*vnumnodes+i];
4859 Bprime[dim*i+2] = vdbasis[2*vnumnodes+i];
4864 xDelete<IssmDouble>(vdbasis);
◆ GetBFSprimevel()
Definition at line 4866 of file StressbalanceAnalysis.cpp.
4889 IssmDouble* vdbasis=xNew<IssmDouble>(dim*vnumnodes);
4894 for(
int i=0;i<vnumnodes;i++){
4895 Bprime[(dim*vnumnodes)*0+dim*i+0] = vdbasis[0*vnumnodes+i];
4896 Bprime[(dim*vnumnodes)*0+dim*i+1] = 0.;
4897 Bprime[(dim*vnumnodes)*1+dim*i+0] = 0.;
4898 Bprime[(dim*vnumnodes)*1+dim*i+1] = vdbasis[1*vnumnodes+i];
4899 Bprime[(dim*vnumnodes)*2+dim*i+0] = vdbasis[1*vnumnodes+i];
4900 Bprime[(dim*vnumnodes)*2+dim*i+1] = vdbasis[0*vnumnodes+i];
4904 for(
int i=0;i<vnumnodes;i++){
4905 Bprime[(dim*vnumnodes)*0+dim*i+0] = vdbasis[0*vnumnodes+i];
4906 Bprime[(dim*vnumnodes)*0+dim*i+1] = 0.;
4907 Bprime[(dim*vnumnodes)*0+dim*i+2] = 0.;
4908 Bprime[(dim*vnumnodes)*1+dim*i+0] = 0.;
4909 Bprime[(dim*vnumnodes)*1+dim*i+1] = vdbasis[1*vnumnodes+i];
4910 Bprime[(dim*vnumnodes)*1+dim*i+2] = 0.;
4911 Bprime[(dim*vnumnodes)*2+dim*i+0] = 0.;
4912 Bprime[(dim*vnumnodes)*2+dim*i+1] = 0.;
4913 Bprime[(dim*vnumnodes)*2+dim*i+2] = vdbasis[2*vnumnodes+i];
4914 Bprime[(dim*vnumnodes)*3+dim*i+0] = vdbasis[1*vnumnodes+i];
4915 Bprime[(dim*vnumnodes)*3+dim*i+1] = vdbasis[0*vnumnodes+i];
4916 Bprime[(dim*vnumnodes)*3+dim*i+2] = 0.;
4917 Bprime[(dim*vnumnodes)*4+dim*i+0] = vdbasis[2*vnumnodes+i];
4918 Bprime[(dim*vnumnodes)*4+dim*i+1] = 0.;
4919 Bprime[(dim*vnumnodes)*4+dim*i+2] = vdbasis[0*vnumnodes+i];
4920 Bprime[(dim*vnumnodes)*5+dim*i+0] = 0.;
4921 Bprime[(dim*vnumnodes)*5+dim*i+1] = vdbasis[2*vnumnodes+i];
4922 Bprime[(dim*vnumnodes)*5+dim*i+2] = vdbasis[1*vnumnodes+i];
4927 xDelete<IssmDouble>(vdbasis);
◆ GetBFSUzawa()
Definition at line 4929 of file StressbalanceAnalysis.cpp.
4935 if(dim==2) pnumnodes=3;
4940 IssmDouble* basis =xNew<IssmDouble>(pnumnodes);
4944 for(
int i=0;i<pnumnodes;i++){
4949 xDelete<IssmDouble>(basis);
◆ GetBFSvel()
Definition at line 4951 of file StressbalanceAnalysis.cpp.
4975 IssmDouble* vdbasis=xNew<IssmDouble>(dim*vnumnodes);
4980 for(
int i=0;i<vnumnodes;i++){
4981 B[(dim*vnumnodes)*0+dim*i+0] = vdbasis[0*vnumnodes+i];
4982 B[(dim*vnumnodes)*0+dim*i+1] = 0.;
4983 B[(dim*vnumnodes)*1+dim*i+0] = 0.;
4984 B[(dim*vnumnodes)*1+dim*i+1] = vdbasis[1*vnumnodes+i];
4985 B[(dim*vnumnodes)*2+dim*i+0] = .5*vdbasis[1*vnumnodes+i];
4986 B[(dim*vnumnodes)*2+dim*i+1] = .5*vdbasis[0*vnumnodes+i];
4990 for(
int i=0;i<vnumnodes;i++){
4991 B[(dim*vnumnodes)*0+dim*i+0] = vdbasis[0*vnumnodes+i];
4992 B[(dim*vnumnodes)*0+dim*i+1] = 0.;
4993 B[(dim*vnumnodes)*0+dim*i+2] = 0.;
4994 B[(dim*vnumnodes)*1+dim*i+0] = 0.;
4995 B[(dim*vnumnodes)*1+dim*i+1] = vdbasis[1*vnumnodes+i];
4996 B[(dim*vnumnodes)*1+dim*i+2] = 0.;
4997 B[(dim*vnumnodes)*2+dim*i+0] = 0.;
4998 B[(dim*vnumnodes)*2+dim*i+1] = 0.;
4999 B[(dim*vnumnodes)*2+dim*i+2] = vdbasis[2*vnumnodes+i];
5000 B[(dim*vnumnodes)*3+dim*i+0] = .5*vdbasis[1*vnumnodes+i];
5001 B[(dim*vnumnodes)*3+dim*i+1] = .5*vdbasis[0*vnumnodes+i];
5002 B[(dim*vnumnodes)*3+dim*i+2] = 0.;
5003 B[(dim*vnumnodes)*4+dim*i+0] = .5*vdbasis[2*vnumnodes+i];
5004 B[(dim*vnumnodes)*4+dim*i+1] = 0.;
5005 B[(dim*vnumnodes)*4+dim*i+2] = .5*vdbasis[0*vnumnodes+i];
5006 B[(dim*vnumnodes)*5+dim*i+0] = 0.;
5007 B[(dim*vnumnodes)*5+dim*i+1] = .5*vdbasis[2*vnumnodes+i];
5008 B[(dim*vnumnodes)*5+dim*i+2] = .5*vdbasis[1*vnumnodes+i];
5013 xDelete<IssmDouble>(vdbasis);
◆ GetCFS()
Definition at line 5015 of file StressbalanceAnalysis.cpp.
5024 IssmDouble* basis =xNew<IssmDouble>(lnumnodes);
5028 for(
int i=0;i<lnumnodes;i++){
5029 C[lnumnodes*0+i] = basis[i];
5030 C[lnumnodes*1+i] = basis[i];
5034 xDelete<IssmDouble>(basis);
◆ GetCFSprime()
Definition at line 5036 of file StressbalanceAnalysis.cpp.
5050 int vnumdof = vnumnodes*dim;
5052 IssmDouble* vbasis=xNew<IssmDouble>(vnumnodes);
5057 for(
int i=0;i<vnumnodes;i++){
5058 Cprime[vnumdof*0+3*i+0] = vbasis[i];
5059 Cprime[vnumdof*0+3*i+1] = 0.;
5060 Cprime[vnumdof*0+3*i+2] = 0.;
5062 Cprime[vnumdof*1+3*i+0] = 0.;
5063 Cprime[vnumdof*1+3*i+1] = vbasis[i];
5064 Cprime[vnumdof*1+3*i+2] = 0.;
5066 Cprime[vnumdof*2+3*i+0] = 0.;
5067 Cprime[vnumdof*2+3*i+1] = 0.;
5068 Cprime[vnumdof*2+3*i+2] = vbasis[i];
5072 for(
int i=0;i<vnumnodes;i++){
5073 Cprime[vnumdof*0+2*i+0] = vbasis[i];
5074 Cprime[vnumdof*0+2*i+1] = 0.;
5076 Cprime[vnumdof*1+2*i+0] = 0.;
5077 Cprime[vnumdof*1+2*i+1] = vbasis[i];
5082 xDelete<IssmDouble>(vbasis);
◆ GetSolutionFromInputsFS()
Definition at line 5084 of file StressbalanceAnalysis.cpp.
5100 int vnumdof = vnumnodes*dim;
5101 int pnumdof = pnumnodes*1;
5104 IssmDouble* vvalues = xNew<IssmDouble>(vnumdof);
5105 IssmDouble* pvalues = xNew<IssmDouble>(pnumdof);
5119 for(
int i=0;i<vnumnodes;i++){
5123 vvalues[i*dim+0]=vx;
5124 vvalues[i*dim+1]=vy;
5127 vvalues[i*dim+2]=vz;
5130 for(
int i=0;i<pnumnodes;i++){
5133 pvalues[i]=p/FSreconditioning;
5139 vvalues[(vnumnodes-1)*dim+0]=0.;
5140 vvalues[(vnumnodes-1)*dim+1]=0.;
5141 if(dim==3) vvalues[(vnumnodes-1)*dim+2]=0.;
5150 xDelete<int>(pdoflist);
5151 xDelete<int>(vdoflist);
5152 xDelete<IssmDouble>(pvalues);
5153 xDelete<IssmDouble>(vvalues);
◆ InitializeXTH()
void StressbalanceAnalysis::InitializeXTH |
( |
Elements * |
elements, |
|
|
Parameters * |
parameters |
|
) |
| |
Definition at line 5155 of file StressbalanceAnalysis.cpp.
5166 for(
int i=0;i<elements->
Size();i++){
5178 IssmDouble* epsxx = xNew<IssmDouble>(tnumnodes);
IssmDouble* sigmapxx = xNew<IssmDouble>(tnumnodes);
5179 IssmDouble* epsyy = xNew<IssmDouble>(tnumnodes);
IssmDouble* sigmapyy = xNew<IssmDouble>(tnumnodes);
5180 IssmDouble* epsxy = xNew<IssmDouble>(tnumnodes);
IssmDouble* sigmapxy = xNew<IssmDouble>(tnumnodes);
5185 epszz = xNew<IssmDouble>(tnumnodes); sigmapzz = xNew<IssmDouble>(tnumnodes);
5186 epsxz = xNew<IssmDouble>(tnumnodes); sigmapxz = xNew<IssmDouble>(tnumnodes);
5187 epsyz = xNew<IssmDouble>(tnumnodes); sigmapyz = xNew<IssmDouble>(tnumnodes);
5192 for(
int i=0;i<tnumnodes;i++){
5202 epsxx[i] = dvx[0]; sigmapxx[i] = 2.*viscosity*epsxx[i];
5203 epsyy[i] = dvy[1]; sigmapyy[i] = 2.*viscosity*epsyy[i];
5204 epsxy[i] = 0.5*(dvx[1] + dvy[0]); sigmapxy[i] = 2.*viscosity*epsxy[i];
5206 epszz[i] = dvz[2]; sigmapzz[i] = 2.*viscosity*epszz[i];
5207 epsxz[i] = 0.5*(dvx[2] + dvz[0]); sigmapxz[i] = 2.*viscosity*epsxz[i];
5208 epsyz[i] = 0.5*(dvy[2] + dvz[1]); sigmapyz[i] = 2.*viscosity*epsyz[i];
5224 xDelete<IssmDouble>(xyz_list);
5225 xDelete<IssmDouble>(epsxx); xDelete<IssmDouble>(sigmapxx);
5226 xDelete<IssmDouble>(epsyy); xDelete<IssmDouble>(sigmapyy);
5227 xDelete<IssmDouble>(epszz); xDelete<IssmDouble>(sigmapzz);
5228 xDelete<IssmDouble>(epsxy); xDelete<IssmDouble>(sigmapxy);
5229 xDelete<IssmDouble>(epsxz); xDelete<IssmDouble>(sigmapxz);
5230 xDelete<IssmDouble>(epsyz); xDelete<IssmDouble>(sigmapyz);
◆ InputUpdateFromSolutionFS()
void StressbalanceAnalysis::InputUpdateFromSolutionFS |
( |
IssmDouble * |
solution, |
|
|
Element * |
element |
|
) |
| |
Definition at line 5234 of file StressbalanceAnalysis.cpp.
5247 int vnumdof = vnumnodes*dim;
5248 int pnumdof = pnumnodes*1;
5251 IssmDouble* values = xNew<IssmDouble>(vnumdof+pnumdof);
5252 IssmDouble* vx = xNew<IssmDouble>(vnumnodes);
5253 IssmDouble* vy = xNew<IssmDouble>(vnumnodes);
5254 IssmDouble* vz = xNew<IssmDouble>(vnumnodes);
5255 IssmDouble* vel = xNew<IssmDouble>(vnumnodes);
5256 IssmDouble* pressure = xNew<IssmDouble>(pnumnodes);
5259 int* cs_list = xNew<int>(vnumnodes+pnumnodes);
5260 if(dim==2)
for(i=0;i<vnumnodes;i++) cs_list[i] =
XYEnum;
5261 else for(i=0;i<vnumnodes;i++) cs_list[i] =
XYZEnum;
5262 for(i=0;i<pnumnodes;i++) cs_list[vnumnodes+i] =
PressureEnum;
5269 for(i=0;i<vnumdof;i++) values[i] =solution[vdoflist[i]];
5270 for(i=0;i<pnumdof;i++) values[vnumdof+i]=solution[pdoflist[i]];
5276 for(i=0;i<vnumnodes;i++){
5277 vx[i] = values[i*dim+0];
5278 vy[i] = values[i*dim+1];
5279 if(xIsNan<IssmDouble>(vx[i]))
_error_(
"NaN found in solution vector");
5280 if(xIsInf<IssmDouble>(vx[i]))
_error_(
"Inf found in solution vector");
5281 if(xIsNan<IssmDouble>(vy[i]))
_error_(
"NaN found in solution vector");
5282 if(xIsInf<IssmDouble>(vy[i]))
_error_(
"Inf found in solution vector");
5285 vz[i] = values[i*dim+2];
5286 if(xIsNan<IssmDouble>(vz[i]))
_error_(
"NaN found in solution vector");
5287 if(xIsInf<IssmDouble>(vz[i]))
_error_(
"Inf found in solution vector");
5290 for(i=0;i<pnumnodes;i++){
5291 pressure[i] = values[vnumdof+i];
5292 if(xIsNan<IssmDouble>(pressure[i]))
_error_(
"NaN found in solution vector");
5293 if(xIsInf<IssmDouble>(pressure[i]))
_error_(
"Inf found in solution vector");
5297 for(i=0;i<pnumnodes;i++) pressure[i] = pressure[i]*FSreconditioning;
5298 if(dim==3)
for(i=0;i<vnumnodes;i++) vel[i] = sqrt(vx[i]*vx[i] + vy[i]*vy[i] + vz[i]*vz[i]);
5299 else for(i=0;i<vnumnodes;i++) vel[i] = sqrt(vx[i]*vx[i] + vy[i]*vy[i]);
5312 xDelete<IssmDouble>(pressure);
5313 xDelete<IssmDouble>(vel);
5314 xDelete<IssmDouble>(vz);
5315 xDelete<IssmDouble>(vy);
5316 xDelete<IssmDouble>(vx);
5317 xDelete<IssmDouble>(values);
5318 xDelete<int>(vdoflist);
5319 xDelete<int>(pdoflist);
5320 xDelete<int>(cs_list);
◆ InputUpdateFromSolutionFSXTH_d()
void StressbalanceAnalysis::InputUpdateFromSolutionFSXTH_d |
( |
Elements * |
elements, |
|
|
Parameters * |
parameters |
|
) |
| |
Definition at line 5322 of file StressbalanceAnalysis.cpp.
5326 IssmDouble epsxx,epsyy,epszz,epsxy,epsxz,epsyz,D_scalar;
5327 IssmDouble epsxx_old,epsyy_old,epszz_old,epsxy_old,epsxz_old,epsyz_old;
5328 IssmDouble sigmapxx,sigmapyy,sigmapzz,sigmapxy,sigmapxz,sigmapyz;
5335 if(dim==2) tausize = 3;
5338 for(
int i=0;i<elements->
Size();i++){
5354 IssmDouble* tbasis = xNew<IssmDouble>(tnumnodes);
5355 IssmDouble* Ke = xNewZeroInit<IssmDouble>(tnumnodes*tnumnodes);
5356 IssmDouble* pe_xx = xNewZeroInit<IssmDouble>(tnumnodes);
5357 IssmDouble* pe_yy = xNewZeroInit<IssmDouble>(tnumnodes);
5358 IssmDouble* pe_xy = xNewZeroInit<IssmDouble>(tnumnodes);
5361 pe_zz = xNewZeroInit<IssmDouble>(tnumnodes);
5362 pe_xz = xNewZeroInit<IssmDouble>(tnumnodes);
5363 pe_yz = xNewZeroInit<IssmDouble>(tnumnodes);
5381 Input2* sigmapzz_input=NULL;
Input2* sigmapxz_input=NULL;
Input2* sigmapyz_input=NULL;
5389 for(
int ig=gauss->
begin();ig<gauss->end();ig++){
5422 epsxy = 0.5*(dvx[1] + dvy[0]);
5425 epsxz = 0.5*(dvx[2] + dvz[0]);
5426 epsyz = 0.5*(dvy[2] + dvz[1]);
5433 coef1 = B*pow(1./sqrt(2.),(1.-n)/n);
5437 (r*epsxx + sigmapxx)*(r*epsxx + sigmapxx)
5438 + (r*epsyy + sigmapyy)*(r*epsyy + sigmapyy)
5439 + 2*(r*epsxy + sigmapxy)*(r*epsxy + sigmapxy)
5444 (r*epsxx + sigmapxx)*(r*epsxx + sigmapxx)
5445 + (r*epsyy + sigmapyy)*(r*epsyy + sigmapyy)
5446 + (r*epszz + sigmapzz)*(r*epszz + sigmapzz)
5447 + 2*(r*epsxy + sigmapxy)*(r*epsxy + sigmapxy)
5448 + 2*(r*epsxz + sigmapxz)*(r*epsxz + sigmapxz)
5449 + 2*(r*epsyz + sigmapyz)*(r*epsyz + sigmapyz)
5454 dnorm = sqrt( epsxx_old*epsxx_old + epsyy_old*epsyy_old + 2.*epsxy_old*epsxy_old );
5457 dnorm = sqrt( epsxx_old*epsxx_old + epsyy_old*epsyy_old + epszz_old*epszz_old
5458 +2.*(epsxy_old*epsxy_old + epsxz_old*epsxz_old + epsyz_old*epsyz_old));
5461 if(dnorm==0.) dnorm=1.;
5464 _assert_(!xIsNan<IssmDouble>(dnorm));
5467 D_scalar=(coef1*pow(dnorm,(1.-n)/n)+r)*gauss->
weight*Jdet;
5470 tbasis,1,tnumnodes,0,
5474 for(
int ii=0;ii<tnumnodes;ii++) pe_xx[ii] += (r*epsxx+sigmapxx)*tbasis[ii]*gauss->
weight*Jdet;
5475 for(
int ii=0;ii<tnumnodes;ii++) pe_yy[ii] += (r*epsyy+sigmapyy)*tbasis[ii]*gauss->
weight*Jdet;
5476 for(
int ii=0;ii<tnumnodes;ii++) pe_xy[ii] += (r*epsxy+sigmapxy)*tbasis[ii]*gauss->
weight*Jdet;
5478 for(
int ii=0;ii<tnumnodes;ii++) pe_zz[ii] += (r*epszz+sigmapzz)*tbasis[ii]*gauss->
weight*Jdet;
5479 for(
int ii=0;ii<tnumnodes;ii++) pe_xz[ii] += (r*epsxz+sigmapxz)*tbasis[ii]*gauss->
weight*Jdet;
5480 for(
int ii=0;ii<tnumnodes;ii++) pe_yz[ii] += (r*epsyz+sigmapyz)*tbasis[ii]*gauss->
weight*Jdet;
5485 IssmDouble* d_xx = xNew<IssmDouble>(tnumnodes);
5486 IssmDouble* d_yy = xNew<IssmDouble>(tnumnodes);
5487 IssmDouble* d_xy = xNew<IssmDouble>(tnumnodes);
5496 for(
int i=0;i<3;i++) _assert_(!xIsNan<IssmDouble>(d_xx[i]));
5497 for(
int i=0;i<3;i++) _assert_(!xIsNan<IssmDouble>(d_yy[i]));
5498 for(
int i=0;i<3;i++) _assert_(!xIsNan<IssmDouble>(d_xx[i]));
5505 d_zz = xNew<IssmDouble>(tnumnodes);
5506 d_xz = xNew<IssmDouble>(tnumnodes);
5507 d_yz = xNew<IssmDouble>(tnumnodes);
5524 xDelete<IssmDouble>(xyz_list);
5525 xDelete<IssmDouble>(tbasis);
5526 xDelete<IssmDouble>(Ke);
5527 xDelete<IssmDouble>(pe_xx); xDelete<IssmDouble>(d_xx);
5528 xDelete<IssmDouble>(pe_yy); xDelete<IssmDouble>(d_yy);
5529 xDelete<IssmDouble>(pe_zz); xDelete<IssmDouble>(d_zz);
5530 xDelete<IssmDouble>(pe_xy); xDelete<IssmDouble>(d_xy);
5531 xDelete<IssmDouble>(pe_xz); xDelete<IssmDouble>(d_xz);
5532 xDelete<IssmDouble>(pe_yz); xDelete<IssmDouble>(d_yz);
◆ InputUpdateFromSolutionFSXTH_tau()
void StressbalanceAnalysis::InputUpdateFromSolutionFSXTH_tau |
( |
Elements * |
elements, |
|
|
Parameters * |
parameters |
|
) |
| |
Definition at line 5535 of file StressbalanceAnalysis.cpp.
5539 IssmDouble epsxx,epsyy,epszz,epsxy,epsxz,epsyz,D_scalar;
5541 IssmDouble sigmapxx,sigmapyy,sigmapzz,sigmapxy,sigmapxz,sigmapyz;
5548 if(dim==2) tausize = 3;
5551 for(
int i=0;i<elements->
Size();i++){
5565 Input2* sigmapzz_input=NULL;
Input2* sigmapxz_input=NULL;
Input2* sigmapyz_input=NULL;
5587 IssmDouble* tau_xx = xNew<IssmDouble>(tnumnodes);
5588 IssmDouble* tau_yy = xNew<IssmDouble>(tnumnodes);
5589 IssmDouble* tau_xy = xNew<IssmDouble>(tnumnodes);
5594 tau_zz = xNew<IssmDouble>(tnumnodes);
5595 tau_xz = xNew<IssmDouble>(tnumnodes);
5596 tau_yz = xNew<IssmDouble>(tnumnodes);
5599 for(
int ig=0;ig<tnumnodes;ig++){
5610 epsxy = 0.5*(dvx[1] + dvy[0]);
5613 epsxz = 0.5*(dvx[2] + dvz[0]);
5614 epsyz = 0.5*(dvy[2] + dvz[1]);
5638 tau_xx[ig] = sigmapxx + r*(epsxx - d_xx);
5639 tau_yy[ig] = sigmapyy + r*(epsyy - d_yy);
5640 tau_xy[ig] = sigmapxy + r*(epsxy - d_xy);
5642 tau_zz[ig] = sigmapzz + r*(epszz - d_zz);
5643 tau_xz[ig] = sigmapxz + r*(epsxz - d_xz);
5644 tau_yz[ig] = sigmapyz + r*(epsyz - d_yz);
5660 xDelete<IssmDouble>(xyz_list);
5661 xDelete<IssmDouble>(tau_xx);
5662 xDelete<IssmDouble>(tau_yy);
5663 xDelete<IssmDouble>(tau_zz);
5664 xDelete<IssmDouble>(tau_xy);
5665 xDelete<IssmDouble>(tau_xz);
5666 xDelete<IssmDouble>(tau_yz);
◆ CreateKMatrixCouplingHOFS()
Definition at line 5671 of file StressbalanceAnalysis.cpp.
5677 int numnodes = 3*6+1;
5679 int numdofs = 4*6 + 3;
5680 int numdoftotal = (2+4)*6+ 3;
5684 int* cs_list = xNew<int>(6*3+1);
5685 int* cs_list2 = xNew<int>(6*2+1);
5686 Node **node_list = xNew<Node*>(6*3+1);
5693 node_list[i+6] = element->
GetNode(i);
5698 node_list[i] = element->
GetNode(i);
5699 node_list[i+2*6+1] = element->
GetNode(i+6*1);
5709 delete Ke1;
delete Ke2;
5713 int indices[3]={18,19,20};
5720 for(i=0;i<numdofs;i++)
for(j=0;j<6;j++){
5721 Ke->
values[(i+numdofp)*numdoftotal+2*j+0]+=Ke1->
values[i*numdofs+3*j+0];
5722 Ke->
values[(i+numdofp)*numdoftotal+2*j+1]+=Ke1->
values[i*numdofs+3*j+1];
5724 for(i=0;i<numdofp;i++)
for(j=0;j<6;j++){
5725 Ke->
values[i*numdoftotal+numdofp+3*j+0]+=Ke2->
values[i*numdofp+2*j+0];
5726 Ke->
values[i*numdoftotal+numdofp+3*j+1]+=Ke2->
values[i*numdofp+2*j+1];
5733 xDelete<int>(cs_list);
5734 xDelete<int>(cs_list2);
5735 xDelete<Node*>(node_list);
◆ CreateKMatrixCouplingSSAFS()
◆ CreateKMatrixCouplingSSAFSFriction()
Definition at line 5752 of file StressbalanceAnalysis.cpp.
5758 const int numdofs = (6+1)*3 + 6*1;
5759 const int numdofm = 6 *2;
5760 const int numdof2d = 3 *3;
5761 const int numdof2dm = 3 *2;
5762 const int numdoftot = 6*2 + (6+1)*3 +6;
5765 int i,j,approximation;
5767 IssmDouble FSreconditioning,viscosity,alpha2_gauss,Jdet2d;
5775 IssmDouble Ke_drag[numdof2dm][numdofs]={0.0};
5776 IssmDouble Ke_drag2[numdof2d][numdof2dm]={0.0};
5786 int numnodes = 2*vnumnodes-1+pnumnodes;
5789 int* cs_list = xNew<int>(2*vnumnodes-1+pnumnodes);
5790 Node **node_list = xNew<Node*>(2*vnumnodes-1+pnumnodes);
5791 for(i=0;i<vnumnodes-1;i++){
5792 node_list[i] = element->
GetNode(i);
5795 for(i=0;i<vnumnodes;i++){
5796 node_list[i+vnumnodes-1] = element->
GetNode(i);
5797 cs_list[i+vnumnodes-1] =
XYZEnum;
5799 for(i=0;i<pnumnodes;i++){
5800 node_list[2*vnumnodes-1+i] = element->
GetNode(vnumnodes+i);
5807 delete Ke1;
delete Ke2;
5822 for(
int ig=gauss->
begin();ig<gauss->end();ig++){
5827 this->
GetLSSAFS(&LSSAFS[0][0], element,gauss);
5828 this->
GetLprimeSSAFS(&LprimeSSAFS[0][0], element,xyz_list, gauss);
5829 this->
GetLFSSSA(&LFSSSA[0][0],element, gauss);
5830 this->
GetLprimeFSSSA(&LprimeFSSSA[0][0], element,xyz_list, gauss);
5834 element->
NormalBase(&bed_normal[0],xyz_list_tria);
5835 friction->
GetAlpha2(&alpha2_gauss,gauss);
5837 DLSSAFS[0][0]=alpha2_gauss*gauss->
weight*Jdet2d;
5838 DLSSAFS[1][1]=alpha2_gauss*gauss->
weight*Jdet2d;
5839 DLSSAFS[2][2]=-alpha2_gauss*gauss->
weight*Jdet2d*bed_normal[0]*bed_normal[2];
5840 DLSSAFS[3][3]=-alpha2_gauss*gauss->
weight*Jdet2d*bed_normal[1]*bed_normal[2];
5841 DLSSAFS[4][4]=-2*viscosity*gauss->
weight*Jdet2d*bed_normal[0];
5842 DLSSAFS[5][5]=-2*viscosity*gauss->
weight*Jdet2d*bed_normal[1];
5843 DLSSAFS[6][6]=FSreconditioning*gauss->
weight*Jdet2d*bed_normal[0];
5844 DLSSAFS[7][7]=FSreconditioning*gauss->
weight*Jdet2d*bed_normal[1];
5846 DLFSSSA[0][0]=alpha2_gauss*gauss->
weight*Jdet2d;
5847 DLFSSSA[1][1]=alpha2_gauss*gauss->
weight*Jdet2d;
5848 DLFSSSA[2][2]=-alpha2_gauss*gauss->
weight*Jdet2d*bed_normal[0]*bed_normal[2];
5849 DLFSSSA[3][3]=-alpha2_gauss*gauss->
weight*Jdet2d*bed_normal[1]*bed_normal[2];
5852 &DLSSAFS[0][0],8,8,0,
5853 &LprimeSSAFS[0][0],8,numdofs,0,
5857 &DLFSSSA[0][0],4,4,0,
5858 &LprimeFSSSA[0][0],4,numdof2dm,0,
5862 for(i=0;i<numdof2dm;i++)
for(j=0;j<numdofs;j++) Ke->
values[i*numdoftot+j+numdofm]+=Ke_drag[i][j];
5863 for(i=0;i<numdof2d;i++)
for(j=0;j<numdof2dm;j++) Ke->
values[(i+numdofm)*numdoftot+j]+=Ke_drag2[i][j];
5869 xDelete<int>(cs_list);
5870 xDelete<Node*>(node_list);
5871 xDelete<IssmDouble>(xyz_list);
5872 xDelete<IssmDouble>(xyz_list_tria);
◆ CreateKMatrixCouplingSSAFSViscous()
Definition at line 5877 of file StressbalanceAnalysis.cpp.
5883 const int numdofm = 2 *3;
5884 const int numdofs = 4 *6+ 3;
5885 const int numdoftotal = 2 *numdofm+numdofs;
5890 IssmDouble Jdet,viscosity,FSreconditioning,D_scalar;
5907 int numnodes = 2*vnumnodes-1+pnumnodes;
5910 int* cs_list = xNew<int>(2*vnumnodes-1+pnumnodes);
5911 Node **node_list = xNew<Node*>(2*vnumnodes-1+pnumnodes);
5912 for(i=0;i<vnumnodes-1;i++){
5913 node_list[i] = pentabase->
GetNode(i);
5916 for(i=0;i<vnumnodes;i++){
5917 node_list[i+vnumnodes-1] = element->
GetNode(i);
5918 cs_list[i+vnumnodes-1] =
XYZEnum;
5920 for(i=0;i<pnumnodes;i++){
5921 node_list[2*vnumnodes-1+i] = element->
GetNode(vnumnodes+i);
5929 delete Ke1;
delete Ke2;
5941 for(
int ig=gauss->
begin();ig<gauss->end();ig++){
5947 this->
GetBSSAFS(&B[0][0],element,xyz_list, gauss);
5949 this->
GetBSSAFSTria(&B2[0][0], basaltria,xyz_list, gauss_tria);
5954 D_scalar=2*viscosity*gauss->
weight*Jdet;
5955 for (i=0;i<3;i++) D[i][i]=D_scalar;
5956 D[3][3]=-gauss->
weight*Jdet*FSreconditioning;
5957 for (i=0;i<3;i++)
D2[i][i]=D_scalar;
5961 &Bprime[0][0],4,numdofm,0,
5966 &Bprime2[0][0],3,numdofs,0,
5970 for(i=0;i<numdofs;i++)
for(j=0;j<numdofm;j++) Ke->
values[(i+2*numdofm)*numdoftotal+j]+=Ke_gg[i][j];
5971 for(i=0;i<numdofm;i++)
for(j=0;j<numdofs;j++) Ke->
values[i*numdoftotal+(j+2*numdofm)]+=Ke_gg2[i][j];
5977 xDelete<int>(cs_list);
5978 xDelete<Node*>(node_list);
5979 xDelete<IssmDouble>(xyz_list);
5980 delete basaltria->
material;
delete basaltria;
◆ CreateKMatrixCouplingSSAHO()
◆ CreateKMatrixCouplingSSAHOFriction()
Definition at line 5998 of file StressbalanceAnalysis.cpp.
6007 int numdof = 2*numnodes;
6008 int numdoftotal = 4*numnodes;
6014 IssmDouble* L = xNewZeroInit<IssmDouble>(2*numdof);
6017 IssmDouble* Ke_gg = xNewZeroInit<IssmDouble>(numdof*numdof);
6018 Node **node_list = xNew<Node*>(2*numnodes);
6019 int* cs_list = xNew<int>(2*numnodes);
6025 delete Ke1;
delete Ke2;
6028 for(i=0;i<numnodes;i++){
6029 node_list[i+0*numnodes] = element->
GetNode(i);
6030 node_list[i+1*numnodes] = element->
GetNode(i);
6031 cs_list[i+0*numnodes] =
XYEnum;
6032 cs_list[i+1*numnodes] =
XYEnum;
6043 for(
int ig=gauss->
begin();ig<gauss->end();ig++){
6052 DL_scalar=alpha2*gauss->
weight*Jdet2d;
6053 for (i=0;i<2;i++) DL[i][i]=DL_scalar;
6062 for(i=0;i<numdof;i++)
for(j=0;j<numdof;j++) Ke->
values[i*numdoftotal+(numdof+j)]+=Ke_gg[i*numdof+j];
6063 for(i=0;i<numdof;i++)
for(j=0;j<numdof;j++) Ke->
values[(i+numdof)*numdoftotal+j]+=Ke_gg[i*numdof+j];
6071 xDelete<int>(cs_list);
6072 xDelete<Node*>(node_list);
6073 xDelete<IssmDouble>(xyz_list_tria);
6074 xDelete<IssmDouble>(Ke_gg);
6075 xDelete<IssmDouble>(L);
◆ CreateKMatrixCouplingSSAHOViscous()
Definition at line 6078 of file StressbalanceAnalysis.cpp.
6085 int numdofm = 1 *numnodes;
6086 int numdofp = 2 *numnodes;
6087 int numdoftotal = 2 *2 *numnodes;
6094 IssmDouble* Bprime = xNew<IssmDouble>(3*numdofm);
6097 IssmDouble* Ke_gg = xNewZeroInit<IssmDouble>(numdofp*numdofm);
6098 Node **node_list = xNew<Node*>(2*numnodes);
6099 int* cs_list= xNew<int>(2*numnodes);
6106 for(i=0;i<numnodes;i++){
6107 node_list[i+0*numnodes] = pentabase->
GetNode(i);
6108 node_list[i+1*numnodes] = element ->
GetNode(i);
6109 cs_list[i+0*numnodes] =
XYEnum;
6110 cs_list[i+1*numnodes] =
XYEnum;
6117 delete Ke1;
delete Ke2;
6127 for(
int ig=gauss->
begin();ig<gauss->end();ig++){
6133 this->
GetBSSAHO(B, element,xyz_list, gauss);
6134 this->
GetBSSAprime(Bprime,basaltria,2,xyz_list, gauss_tria);
6137 D_scalar=2*viscosity*gauss->
weight*Jdet;
6138 for (i=0;i<3;i++) D[i][i]=D_scalar;
6145 for(i=0;i<numdofp;i++)
for(j=0;j<numdofm;j++) Ke->
values[(i+2*numdofm)*numdoftotal+j]+=Ke_gg[i*numdofm+j];
6146 for(i=0;i<numdofm;i++)
for(j=0;j<numdofp;j++) Ke->
values[i*numdoftotal+(j+2*numdofm)]+=Ke_gg[j*numdofm+i];
6156 xDelete<IssmDouble>(B);
6157 xDelete<IssmDouble>(Bprime);
6158 xDelete<IssmDouble>(Ke_gg);
6159 xDelete<IssmDouble>(xyz_list);
6160 xDelete<Node*>(node_list);
6161 xDelete<int>(cs_list);
◆ CreateKMatrixHOFS()
◆ CreateKMatrixSSAFS()
◆ CreateKMatrixSSAHO()
◆ CreateKMatrixSSA3d()
◆ CreateKMatrixSSA3dFriction()
◆ CreateKMatrixSSA3dViscous()
Definition at line 6247 of file StressbalanceAnalysis.cpp.
6253 const int numdof2d=2*3;
6256 int i,j,approximation;
6285 for(
int ig=gauss->
begin();ig<gauss->end();ig++){
6291 this->
GetBSSA(&B[0][0],basaltria,2,xyz_list, gauss_tria);
6292 this->
GetBSSAprime(&Bprime[0][0],basaltria,2,xyz_list, gauss_tria);
6300 else _error_(
"approximation " << approximation <<
" (" <<
EnumToStringx(approximation) <<
") not supported yet");
6302 D_scalar=2*viscosity*gauss->
weight*Jdet;
6303 for (i=0;i<3;i++) D[i][i]=D_scalar;
6307 &Bprime[0][0],3,numdof2d,0,
6311 for(i=0;i<numdof2d;i++)
for(j=0;j<numdof2d;j++) Ke->
values[i*numdof2d+j]+=Ke_gg[i][j];
6317 xDelete<IssmDouble>(xyz_list);
◆ CreatePVectorSSAFS()
◆ CreatePVectorSSAHO()
◆ CreatePVectorCouplingSSAFS()
◆ CreatePVectorCouplingSSAFSFriction()
Definition at line 6507 of file StressbalanceAnalysis.cpp.
6513 int i,j,approximation;
6531 int* cs_list = xNew<int>(vnumnodes+pnumnodes);
6532 Node **node_list = xNew<Node*>(vnumnodes+pnumnodes);
6533 for(i=0;i<vnumnodes;i++){
6535 node_list[i] = element->
GetNode(i);
6537 for(i=0;i<pnumnodes;i++){
6539 node_list[vnumnodes+i] = element->
GetNode(vnumnodes+i);
6557 for(
int ig=gauss->
begin();ig<gauss->end();ig++){
6567 element->
NormalBase(&bed_normal[0],xyz_list_tria);
6569 friction->
GetAlpha2(&alpha2_gauss,gauss);
6572 pe->
values[i*3+0]+=Jdet2d*gauss->
weight*(alpha2_gauss*w*bed_normal[0]*bed_normal[2]+2*viscosity*dw[2]*bed_normal[0])*basis[i];
6573 pe->
values[i*3+1]+=Jdet2d*gauss->
weight*(alpha2_gauss*w*bed_normal[1]*bed_normal[2]+2*viscosity*dw[2]*bed_normal[1])*basis[i];
6574 pe->
values[i*3+2]+=Jdet2d*gauss->
weight*2*viscosity*(dw[0]*bed_normal[0]+dw[1]*bed_normal[1]+dw[2]*bed_normal[2])*basis[i];
6582 xDelete<int>(cs_list);
6583 xDelete<IssmDouble>(xyz_list);
6584 xDelete<IssmDouble>(xyz_list_tria);
6585 xDelete<Node*>(node_list);
◆ CreatePVectorCouplingSSAFSViscous()
Definition at line 6590 of file StressbalanceAnalysis.cpp.
6596 int i,approximation;
6610 int* cs_list = xNew<int>(vnumnodes+pnumnodes);
6611 Node **node_list = xNew<Node*>(vnumnodes+pnumnodes);
6612 for(i=0;i<vnumnodes;i++){
6614 node_list[i] = element->
GetNode(i);
6616 for(i=0;i<pnumnodes;i++){
6618 node_list[vnumnodes+i] = element->
GetNode(vnumnodes+i);
6632 for(
int ig=gauss->
begin();ig<gauss->end();ig++){
6643 pe->
values[i*3+0]+=-Jdet*gauss->
weight*viscosity*dw[0]*dbasis[2][i];
6644 pe->
values[i*3+1]+=-Jdet*gauss->
weight*viscosity*dw[1]*dbasis[2][i];
6645 pe->
values[i*3+2]+=-Jdet*gauss->
weight*viscosity*(dw[0]*dbasis[0][i]+dw[1]*dbasis[1][i]+2*dw[2]*dbasis[2][i]);
6646 pe->
values[3*vnumnodes+i]+=Jdet*gauss->
weight*FSreconditioning*dw[2]*basis[i];
6654 xDelete<int>(cs_list);
6655 xDelete<Node*>(node_list);
6656 xDelete<IssmDouble>(xyz_list);
◆ CreatePVectorHOFS()
◆ CreatePVectorCouplingHOFS()
◆ CreatePVectorCouplingHOFSFriction()
Definition at line 6337 of file StressbalanceAnalysis.cpp.
6343 int i,approximation;
6360 int numnodes = vnumnodes+pnumnodes;
6363 int* cs_list = xNew<int>(vnumnodes+pnumnodes);
6364 Node **node_list = xNew<Node*>(vnumnodes+pnumnodes);
6365 for(i=0;i<vnumnodes;i++){
6367 node_list[i] = element->
GetNode(i);
6369 for(i=0;i<pnumnodes;i++){
6371 node_list[vnumnodes+i] = element->
GetNode(vnumnodes+i);
6390 for(
int ig=gauss->
begin();ig<gauss->end();ig++){
6400 element->
NormalBase(&bed_normal[0],xyz_list_tria);
6402 friction->
GetAlpha2(&alpha2_gauss,gauss);
6405 pe->
values[i*3+0]+=Jdet2d*gauss->
weight*(alpha2_gauss*w*bed_normal[0]*bed_normal[2]+2*viscosity*dw[2]*bed_normal[0])*basis[i];
6406 pe->
values[i*3+1]+=Jdet2d*gauss->
weight*(alpha2_gauss*w*bed_normal[1]*bed_normal[2]+2*viscosity*dw[2]*bed_normal[1])*basis[i];
6407 pe->
values[i*3+2]+=Jdet2d*gauss->
weight*2*viscosity*(dw[0]*bed_normal[0]+dw[1]*bed_normal[1]+dw[2]*bed_normal[2])*basis[i];
6415 xDelete<int>(cs_list);
6416 xDelete<Node*>(node_list);
6417 xDelete<IssmDouble>(xyz_list);
6418 xDelete<IssmDouble>(xyz_list_tria);
◆ CreatePVectorCouplingHOFSViscous()
Definition at line 6423 of file StressbalanceAnalysis.cpp.
6429 int i,approximation;
6444 int* cs_list = xNew<int>(vnumnodes+pnumnodes);
6445 Node **node_list = xNew<Node*>(vnumnodes+pnumnodes);
6446 for(i=0;i<vnumnodes;i++){
6448 node_list[i] = element->
GetNode(i);
6450 for(i=0;i<pnumnodes;i++){
6452 node_list[vnumnodes+i] = element->
GetNode(vnumnodes+i);
6466 for(
int ig=gauss->
begin();ig<gauss->end();ig++){
6478 pe->
values[i*3+0]+=-Jdet*gauss->
weight*viscosity*dw[0]*dbasis[2][i];
6479 pe->
values[i*3+1]+=-Jdet*gauss->
weight*viscosity*dw[1]*dbasis[2][i];
6480 pe->
values[i*3+2]+=-Jdet*gauss->
weight*viscosity*(dw[0]*dbasis[0][i]+dw[1]*dbasis[1][i]+2*dw[2]*dbasis[2][i]);
6481 pe->
values[3*vnumnodes+i]+=Jdet*gauss->
weight*FSreconditioning*dw[2]*basis[i];
6489 xDelete<int>(cs_list);
6490 xDelete<Node*>(node_list);
6491 xDelete<IssmDouble>(xyz_list);
◆ GetBprimeSSAFS()
Definition at line 6718 of file StressbalanceAnalysis.cpp.
6738 Bprime[(3*7+6)*0+3*i+0] = 2.*dbasismini[0][i];
6739 Bprime[(3*7+6)*0+3*i+1] = dbasismini[1][i];
6740 Bprime[(3*7+6)*0+3*i+2] = 0.;
6741 Bprime[(3*7+6)*1+3*i+0] = dbasismini[0][i];
6742 Bprime[(3*7+6)*1+3*i+1] = 2.*dbasismini[1][i];
6743 Bprime[(3*7+6)*1+3*i+2] = 0.;
6744 Bprime[(3*7+6)*2+3*i+0] = dbasismini[1][i];
6745 Bprime[(3*7+6)*2+3*i+1] = dbasismini[0][i];
6746 Bprime[(3*7+6)*2+3*i+2] = 0.;
6750 Bprime[(3*7+6)*0+3*(6+i)+0] = 0.;
6751 Bprime[(3*7+6)*0+3*(6+i)+1] = 0.;
6752 Bprime[(3*7+6)*0+3*(6+i)+2] = 0.;
6753 Bprime[(3*7+6)*1+3*(6+i)+0] = 0.;
6754 Bprime[(3*7+6)*1+3*(6+i)+1] = 0.;
6755 Bprime[(3*7+6)*1+3*(6+i)+2] = 0.;
6756 Bprime[(3*7+6)*2+3*(6+i)+0] = 0.;
6757 Bprime[(3*7+6)*2+3*(6+i)+1] = 0.;
6758 Bprime[(3*7+6)*2+3*(6+i)+2] = 0.;
6762 Bprime[(3*7+6)*0+7*3+i] = 0.;
6763 Bprime[(3*7+6)*1+7*3+i] = 0.;
6764 Bprime[(3*7+6)*2+7*3+i] = 0.;
◆ GetBprimeSSAFSTria()
Definition at line 6767 of file StressbalanceAnalysis.cpp.
6784 IssmDouble* dbasis=xNew<IssmDouble>(2*numnodes);
6788 for(
int i=0;i<numnodes;i++){
6789 Bprime[2*numnodes*0+2*i+0] = dbasis[0*numnodes+i];
6790 Bprime[2*numnodes*0+2*i+1] = 0.;
6791 Bprime[2*numnodes*1+2*i+0] = 0.;
6792 Bprime[2*numnodes*1+2*i+1] = dbasis[1*numnodes+i];
6793 Bprime[2*numnodes*2+2*i+0] = dbasis[1*numnodes+i];
6794 Bprime[2*numnodes*2+2*i+1] = dbasis[0*numnodes+i];
6795 Bprime[2*numnodes*3+2*i+0] = dbasis[0*numnodes+i];
6796 Bprime[2*numnodes*3+2*i+1] = dbasis[1*numnodes+i];
6800 xDelete<IssmDouble>(dbasis);
◆ GetBSSAFS()
Definition at line 6802 of file StressbalanceAnalysis.cpp.
6825 B[(3*7+6)*0+3*i+0] = dbasismini[0][i];
6826 B[(3*7+6)*0+3*i+1] = 0.;
6827 B[(3*7+6)*0+3*i+2] = 0.;
6828 B[(3*7+6)*1+3*i+0] = 0.;
6829 B[(3*7+6)*1+3*i+1] = dbasismini[1][i];
6830 B[(3*7+6)*1+3*i+2] = 0.;
6831 B[(3*7+6)*2+3*i+0] = 0.5*dbasismini[1][i];
6832 B[(3*7+6)*2+3*i+1] = 0.5*dbasismini[0][i];
6833 B[(3*7+6)*2+3*i+2] = 0.;
6834 B[(3*7+6)*3+3*i+0] = 0.;
6835 B[(3*7+6)*3+3*i+1] = 0.;
6836 B[(3*7+6)*3+3*i+2] = 0.;
6839 B[(3*7+6)*0+3*(6+i)+0] = 0.;
6840 B[(3*7+6)*0+3*(6+i)+1] = 0.;
6841 B[(3*7+6)*0+3*(6+i)+2] = 0.;
6842 B[(3*7+6)*1+3*(6+i)+0] = 0.;
6843 B[(3*7+6)*1+3*(6+i)+1] = 0.;
6844 B[(3*7+6)*1+3*(6+i)+2] = 0.;
6845 B[(3*7+6)*2+3*(6+i)+0] = 0.;
6846 B[(3*7+6)*2+3*(6+i)+1] = 0.;
6847 B[(3*7+6)*2+3*(6+i)+2] = 0.;
6848 B[(3*7+6)*3+3*(6+i)+0] = 0.;
6849 B[(3*7+6)*3+3*(6+i)+1] = 0.;
6850 B[(3*7+6)*3+3*(6+i)+2] = 0.;
6854 B[(3*7+6)*0+7*3+i] = 0;
6855 B[(3*7+6)*1+7*3+i] = 0;
6856 B[(3*7+6)*2+7*3+i] = 0;
6857 B[(3*7+6)*3+7*3+i] = basis[i];
◆ GetBSSAFSTria()
Definition at line 6860 of file StressbalanceAnalysis.cpp.
6876 IssmDouble* dbasis=xNew<IssmDouble>(2*numnodes);
6880 for(
int i=0;i<numnodes;i++){
6881 B[2*numnodes*0+2*i+0] = dbasis[0*numnodes+i];
6882 B[2*numnodes*0+2*i+1] = 0.;
6883 B[2*numnodes*1+2*i+0] = 0.;
6884 B[2*numnodes*1+2*i+1] = dbasis[1*numnodes+i];
6885 B[2*numnodes*2+2*i+0] = 0.5*dbasis[1*numnodes+i];
6886 B[2*numnodes*2+2*i+1] = 0.5*dbasis[0*numnodes+i];
6890 xDelete<IssmDouble>(dbasis);
◆ GetBSSAHO()
Definition at line 6892 of file StressbalanceAnalysis.cpp.
6905 IssmDouble* dbasis=xNew<IssmDouble>(3*numnodes);
6911 for(
int i=0;i<numnodes;i++){
6912 B[2*numnodes*0+2*i+0] = dbasis[0*numnodes+i];
6913 B[2*numnodes*0+2*i+1] = 0.;
6914 B[2*numnodes*1+2*i+0] = 0.;
6915 B[2*numnodes*1+2*i+1] = dbasis[1*numnodes+i];
6916 B[2*numnodes*2+2*i+0] = .5*dbasis[1*numnodes+i];
6917 B[2*numnodes*2+2*i+1] = .5*dbasis[0*numnodes+i];
6921 xDelete<IssmDouble>(dbasis);
◆ GetLFSSSA()
Definition at line 7062 of file StressbalanceAnalysis.cpp.
7078 GaussPenta* gauss = xDynamicCast<GaussPenta*>(gauss_in);
7086 for(
int i=0;i<3;i++){
7087 LFS[num_dof*3*0+num_dof*i+0] = basis[i];
7088 LFS[num_dof*3*0+num_dof*i+1] = 0.;
7089 LFS[num_dof*3*0+num_dof*i+2] = 0.;
7090 LFS[num_dof*3*1+num_dof*i+0] = 0.;
7091 LFS[num_dof*3*1+num_dof*i+1] = basis[i];
7092 LFS[num_dof*3*1+num_dof*i+2] = 0.;
7093 LFS[num_dof*3*2+num_dof*i+0] = 0.;
7094 LFS[num_dof*3*2+num_dof*i+1] = 0.;
7095 LFS[num_dof*3*2+num_dof*i+2] = basis[i];
7096 LFS[num_dof*3*3+num_dof*i+0] = 0.;
7097 LFS[num_dof*3*3+num_dof*i+1] = 0.;
7098 LFS[num_dof*3*3+num_dof*i+2] = basis[i];
◆ GetLprimeFSSSA()
Definition at line 6923 of file StressbalanceAnalysis.cpp.
6938 GaussPenta* gauss = xDynamicCast<GaussPenta*>(gauss_in);
6946 for(
int i=0;i<3;i++){
6947 LprimeFS[num_dof*3*0+num_dof*i+0] = basis[i];
6948 LprimeFS[num_dof*3*0+num_dof*i+1] = 0.;
6949 LprimeFS[num_dof*3*1+num_dof*i+0] = 0.;
6950 LprimeFS[num_dof*3*1+num_dof*i+1] = basis[i];
6951 LprimeFS[num_dof*3*2+num_dof*i+0] = basis[i];
6952 LprimeFS[num_dof*3*2+num_dof*i+1] = 0.;
6953 LprimeFS[num_dof*3*3+num_dof*i+0] = 0.;
6954 LprimeFS[num_dof*3*3+num_dof*i+1] = basis[i];
◆ GetLprimeSSAFS()
Definition at line 6957 of file StressbalanceAnalysis.cpp.
6972 int num_dof_vel=3*7;
6973 int num_dof_total=3*7+1*6;
6979 GaussPenta* gauss = xDynamicCast<GaussPenta*>(gauss_in);
6989 for(
int i=0;i<3;i++){
6990 LprimeFS[num_dof_total*0+num_dof*i+0] = basis[i];
6991 LprimeFS[num_dof_total*0+num_dof*i+1] = 0.;
6992 LprimeFS[num_dof_total*0+num_dof*i+2] = 0.;
6993 LprimeFS[num_dof_total*1+num_dof*i+0] = 0.;
6994 LprimeFS[num_dof_total*1+num_dof*i+1] = basis[i];
6995 LprimeFS[num_dof_total*1+num_dof*i+2] = 0.;
6996 LprimeFS[num_dof_total*2+num_dof*i+0] = 0.;
6997 LprimeFS[num_dof_total*2+num_dof*i+1] = 0.;
6998 LprimeFS[num_dof_total*2+num_dof*i+2] = basis[i];
6999 LprimeFS[num_dof_total*3+num_dof*i+0] = 0.;
7000 LprimeFS[num_dof_total*3+num_dof*i+1] = 0.;
7001 LprimeFS[num_dof_total*3+num_dof*i+2] = basis[i];
7002 LprimeFS[num_dof_total*4+num_dof*i+0] = 0.;
7003 LprimeFS[num_dof_total*4+num_dof*i+1] = 0.;
7004 LprimeFS[num_dof_total*4+num_dof*i+2] = dbasis[2][i];
7005 LprimeFS[num_dof_total*5+num_dof*i+0] = 0.;
7006 LprimeFS[num_dof_total*5+num_dof*i+1] = 0.;
7007 LprimeFS[num_dof_total*5+num_dof*i+2] = dbasis[2][i];
7008 LprimeFS[num_dof_total*6+num_dof*i+0] = 0.;
7009 LprimeFS[num_dof_total*6+num_dof*i+1] = 0.;
7010 LprimeFS[num_dof_total*6+num_dof*i+2] = 0.;
7011 LprimeFS[num_dof_total*7+num_dof*i+0] = 0.;
7012 LprimeFS[num_dof_total*7+num_dof*i+1] = 0.;
7013 LprimeFS[num_dof_total*7+num_dof*i+2] = 0.;
7015 for(
int i=3;i<7;i++){
7016 LprimeFS[num_dof_total*0+num_dof*i+0] = 0.;
7017 LprimeFS[num_dof_total*0+num_dof*i+1] = 0.;
7018 LprimeFS[num_dof_total*0+num_dof*i+2] = 0.;
7019 LprimeFS[num_dof_total*1+num_dof*i+0] = 0.;
7020 LprimeFS[num_dof_total*1+num_dof*i+1] = 0.;
7021 LprimeFS[num_dof_total*1+num_dof*i+2] = 0.;
7022 LprimeFS[num_dof_total*2+num_dof*i+0] = 0.;
7023 LprimeFS[num_dof_total*2+num_dof*i+1] = 0.;
7024 LprimeFS[num_dof_total*2+num_dof*i+2] = 0.;
7025 LprimeFS[num_dof_total*3+num_dof*i+0] = 0.;
7026 LprimeFS[num_dof_total*3+num_dof*i+1] = 0.;
7027 LprimeFS[num_dof_total*3+num_dof*i+2] = 0.;
7028 LprimeFS[num_dof_total*4+num_dof*i+0] = 0.;
7029 LprimeFS[num_dof_total*4+num_dof*i+1] = 0.;
7030 LprimeFS[num_dof_total*4+num_dof*i+2] = 0.;
7031 LprimeFS[num_dof_total*5+num_dof*i+0] = 0.;
7032 LprimeFS[num_dof_total*5+num_dof*i+1] = 0.;
7033 LprimeFS[num_dof_total*5+num_dof*i+2] = 0.;
7034 LprimeFS[num_dof_total*6+num_dof*i+0] = 0.;
7035 LprimeFS[num_dof_total*6+num_dof*i+1] = 0.;
7036 LprimeFS[num_dof_total*6+num_dof*i+2] = 0.;
7037 LprimeFS[num_dof_total*7+num_dof*i+0] = 0.;
7038 LprimeFS[num_dof_total*7+num_dof*i+1] = 0.;
7039 LprimeFS[num_dof_total*7+num_dof*i+2] = 0.;
7041 for(
int i=0;i<3;i++){
7042 LprimeFS[num_dof_total*0+num_dof_vel+i] = 0.;
7043 LprimeFS[num_dof_total*1+num_dof_vel+i] = 0.;
7044 LprimeFS[num_dof_total*2+num_dof_vel+i] = 0.;
7045 LprimeFS[num_dof_total*3+num_dof_vel+i] = 0.;
7046 LprimeFS[num_dof_total*4+num_dof_vel+i] = 0.;
7047 LprimeFS[num_dof_total*5+num_dof_vel+i] = 0.;
7048 LprimeFS[num_dof_total*6+num_dof_vel+i] = basis[i];
7049 LprimeFS[num_dof_total*7+num_dof_vel+i] = basis[i];
7051 for(
int i=3;i<6;i++){
7052 LprimeFS[num_dof_total*0+num_dof_vel+i] = 0.;
7053 LprimeFS[num_dof_total*1+num_dof_vel+i] = 0.;
7054 LprimeFS[num_dof_total*2+num_dof_vel+i] = 0.;
7055 LprimeFS[num_dof_total*3+num_dof_vel+i] = 0.;
7056 LprimeFS[num_dof_total*4+num_dof_vel+i] = 0.;
7057 LprimeFS[num_dof_total*5+num_dof_vel+i] = 0.;
7058 LprimeFS[num_dof_total*6+num_dof_vel+i] = 0.;
7059 LprimeFS[num_dof_total*7+num_dof_vel+i] = 0.;
◆ GetLSSAFS()
Definition at line 7101 of file StressbalanceAnalysis.cpp.
7122 GaussPenta* gauss = xDynamicCast<GaussPenta*>(gauss_in);
7130 for(
int i=0;i<3;i++){
7131 LFS[num_dof*3*0+num_dof*i+0] = basis[i];
7132 LFS[num_dof*3*0+num_dof*i+1] = 0;
7133 LFS[num_dof*3*1+num_dof*i+0] = 0;
7134 LFS[num_dof*3*1+num_dof*i+1] = basis[i];
7135 LFS[num_dof*3*2+num_dof*i+0] = basis[i];
7136 LFS[num_dof*3*2+num_dof*i+1] = 0;
7137 LFS[num_dof*3*3+num_dof*i+0] = 0;
7138 LFS[num_dof*3*3+num_dof*i+1] = basis[i];
7139 LFS[num_dof*3*4+num_dof*i+0] = basis[i];
7140 LFS[num_dof*3*4+num_dof*i+1] = 0;
7141 LFS[num_dof*3*5+num_dof*i+0] = 0;
7142 LFS[num_dof*3*5+num_dof*i+1] = basis[i];
7143 LFS[num_dof*3*6+num_dof*i+0] = basis[i];
7144 LFS[num_dof*3*6+num_dof*i+1] = 0;
7145 LFS[num_dof*3*7+num_dof*i+0] = 0;
7146 LFS[num_dof*3*7+num_dof*i+1] = basis[i];
◆ InputUpdateFromSolutionHOFS()
void StressbalanceAnalysis::InputUpdateFromSolutionHOFS |
( |
IssmDouble * |
solution, |
|
|
Element * |
element |
|
) |
| |
Definition at line 7149 of file StressbalanceAnalysis.cpp.
7153 int* doflistHO = NULL;
7154 int* doflistFSv = NULL;
7155 int* doflistFSp = NULL;
7163 int numdofFSv = 6*3;
7170 IssmDouble* HOvalues = xNew<IssmDouble>(numdofHO);
7171 IssmDouble* FSvalues = xNew<IssmDouble>(numdofFSv+numdofFSp);
7175 IssmDouble* vzHO = xNew<IssmDouble>(numnodes);
7176 IssmDouble* vzFS = xNew<IssmDouble>(numnodes);
7177 IssmDouble* vel = xNew<IssmDouble>(numnodes);
7178 IssmDouble* pressure = xNew<IssmDouble>(numnodes);
7181 int* cs_list = xNew<int>(2*numnodes);
7182 for(i=0;i<numnodes;i++) cs_list[i] =
XYZEnum;
7183 for(i=0;i<numnodes;i++) cs_list[numnodes+i] =
PressureEnum;
7187 for(i=0;i<numdofHO ;i++) HOvalues[i]=solution[doflistHO[i]];
7188 for(i=0;i<numdofFSv;i++) FSvalues[i]=solution[doflistFSv[i]];
7189 for(i=0;i<numdofFSp;i++) FSvalues[numdofFSv+i]=solution[doflistFSp[i]];
7196 for(i=0;i<numnodes;i++){
7197 vx[i] = FSvalues[i*3+0]+HOvalues[i*2+0];
7198 vy[i] = FSvalues[i*3+1]+HOvalues[i*2+1];
7199 vzFS[i] = FSvalues[i*3+2];
7200 pressure[i] = FSvalues[numnodes*3+i]*FSreconditioning;
7203 if(xIsNan<IssmDouble>(vx[i]))
_error_(
"NaN found in solution vector");
7204 if(xIsInf<IssmDouble>(vx[i]))
_error_(
"Inf found in solution vector");
7205 if(xIsNan<IssmDouble>(vy[i]))
_error_(
"NaN found in solution vector");
7206 if(xIsInf<IssmDouble>(vy[i]))
_error_(
"Inf found in solution vector");
7207 if(xIsNan<IssmDouble>(vzFS[i]))
_error_(
"NaN found in solution vector");
7208 if(xIsInf<IssmDouble>(vzFS[i]))
_error_(
"Inf found in solution vector");
7209 if(xIsNan<IssmDouble>(pressure[i]))
_error_(
"NaN found in solution vector");
7210 if(xIsInf<IssmDouble>(pressure[i]))
_error_(
"Inf found in solution vector");
7215 for(i=0;i<numnodes;i++){
7216 vz[i] = vzHO[i]+vzFS[i];
7217 vel[i]= sqrt(vx[i]*vx[i] + vy[i]*vy[i] + vz[i]*vz[i]);
7229 xDelete<IssmDouble>(pressure);
7230 xDelete<IssmDouble>(vel);
7231 xDelete<IssmDouble>(vz);
7232 xDelete<IssmDouble>(vzHO);
7233 xDelete<IssmDouble>(vzFS);
7234 xDelete<IssmDouble>(vy);
7235 xDelete<IssmDouble>(vx);
7236 xDelete<IssmDouble>(FSvalues);
7237 xDelete<IssmDouble>(HOvalues);
7238 xDelete<int>(doflistFSp);
7239 xDelete<int>(doflistFSv);
7240 xDelete<int>(doflistHO);
7241 xDelete<int>(cs_list);
◆ InputUpdateFromSolutionSSAFS()
void StressbalanceAnalysis::InputUpdateFromSolutionSSAFS |
( |
IssmDouble * |
solution, |
|
|
Element * |
element |
|
) |
| |
Definition at line 7243 of file StressbalanceAnalysis.cpp.
7247 int* doflistSSA = NULL;
7248 int* doflistFSv = NULL;
7249 int* doflistFSp = NULL;
7260 int numdof2d = numnodes;
7261 int numdofSSA = 6*2;
7262 int numdofFSv = 6*3;
7269 IssmDouble* SSAvalues = xNew<IssmDouble>(numdofSSA);
7270 IssmDouble* FSvalues = xNew<IssmDouble>(numdofFSv+numdofFSp);
7274 IssmDouble* vzSSA = xNew<IssmDouble>(numnodes);
7275 IssmDouble* vzFS = xNew<IssmDouble>(numnodes);
7276 IssmDouble* vel = xNew<IssmDouble>(numnodes);
7277 IssmDouble* pressure = xNew<IssmDouble>(numnodes);
7280 int* cs_list = xNew<int>(2*numnodes);
7281 for(i=0;i<numnodes;i++) cs_list[i] =
XYZEnum;
7282 for(i=0;i<numnodes;i++) cs_list[numnodes+i] =
PressureEnum;
7286 for(i=0;i<numdof2d;i++){
7287 SSAvalues[i] = solution[doflistSSA[i]];
7288 SSAvalues[i+numdof2d] = solution[doflistSSA[i]];
7290 for(i=0;i<numdofFSv;i++) FSvalues[i]=solution[doflistFSv[i]];
7291 for(i=0;i<numdofFSp;i++) FSvalues[numdofFSv+i]=solution[doflistFSp[i]];
7299 for(i=0;i<numnodes;i++){
7300 vx[i] = FSvalues[i*3+0]+SSAvalues[i*2+0];
7301 vy[i] = FSvalues[i*3+1]+SSAvalues[i*2+1];
7302 vzFS[i] = FSvalues[i*3+2];
7303 pressure[i] = FSvalues[numnodes*3+i]*FSreconditioning;
7306 if(xIsNan<IssmDouble>(vx[i]))
_error_(
"NaN found in solution vector");
7307 if(xIsInf<IssmDouble>(vx[i]))
_error_(
"Inf found in solution vector");
7308 if(xIsNan<IssmDouble>(vy[i]))
_error_(
"NaN found in solution vector");
7309 if(xIsInf<IssmDouble>(vy[i]))
_error_(
"Inf found in solution vector");
7310 if(xIsNan<IssmDouble>(vzFS[i]))
_error_(
"NaN found in solution vector");
7311 if(xIsInf<IssmDouble>(vzFS[i]))
_error_(
"Inf found in solution vector");
7312 if(xIsNan<IssmDouble>(pressure[i]))
_error_(
"NaN found in solution vector");
7313 if(xIsInf<IssmDouble>(pressure[i]))
_error_(
"Inf found in solution vector");
7318 for(i=0;i<numnodes;i++){
7319 vz[i] = vzSSA[i]+vzFS[i];
7320 vel[i]= sqrt(vx[i]*vx[i] + vy[i]*vy[i] + vz[i]*vz[i]);
7332 xDelete<IssmDouble>(pressure);
7333 xDelete<IssmDouble>(vel);
7334 xDelete<IssmDouble>(vz);
7335 xDelete<IssmDouble>(vzSSA);
7336 xDelete<IssmDouble>(vzFS);
7337 xDelete<IssmDouble>(vy);
7338 xDelete<IssmDouble>(vx);
7339 xDelete<IssmDouble>(FSvalues);
7340 xDelete<IssmDouble>(SSAvalues);
7341 xDelete<int>(doflistFSp);
7342 xDelete<int>(doflistFSv);
7343 xDelete<int>(doflistSSA);
7344 xDelete<int>(cs_list);
◆ InputUpdateFromSolutionSSAHO()
void StressbalanceAnalysis::InputUpdateFromSolutionSSAHO |
( |
IssmDouble * |
solution, |
|
|
Element * |
element |
|
) |
| |
Definition at line 7346 of file StressbalanceAnalysis.cpp.
7350 int* SSAdoflist = NULL;
7351 int* HOdoflist = NULL;
7363 int numdof = numnodes*2;
7364 int numdof2d = numnodes;
7369 IssmDouble* HOvalues = xNew<IssmDouble>(numdof);
7370 IssmDouble* SSAvalues = xNew<IssmDouble>(numdof);
7374 IssmDouble* vel = xNew<IssmDouble>(numnodes);
7375 IssmDouble* pressure = xNew<IssmDouble>(numnodes);
7376 IssmDouble* surface = xNew<IssmDouble>(numnodes);
7379 for(i=0;i<numdof2d;i++){
7380 HOvalues[i] = solution[HOdoflist[i]];
7381 SSAvalues[i] = solution[SSAdoflist[i]];
7383 for(i=numdof2d;i<numdof;i++){
7384 HOvalues[i] = solution[HOdoflist[i]];
7385 SSAvalues[i] = SSAvalues[i-numdof2d];
7393 for(i=0;i<numnodes;i++){
7394 vx[i]=SSAvalues[i*2+0]+HOvalues[i*2+0];
7395 vy[i]=SSAvalues[i*2+1]+HOvalues[i*2+1];
7398 if(xIsNan<IssmDouble>(vx[i]))
_error_(
"NaN found in solution vector");
7399 if(xIsInf<IssmDouble>(vx[i]))
_error_(
"Inf found in solution vector");
7400 if(xIsNan<IssmDouble>(vy[i]))
_error_(
"NaN found in solution vector");
7401 if(xIsInf<IssmDouble>(vy[i]))
_error_(
"Inf found in solution vector");
7406 for(i=0;i<numnodes;i++) vel[i]=sqrt(vx[i]*vx[i] + vy[i]*vy[i] + vz[i]*vz[i]);
7414 for(i=0;i<numnodes;i++) pressure[i]=rho_ice*g*(surface[i]-xyz_list[i*3+2]);
7423 xDelete<IssmDouble>(surface);
7424 xDelete<IssmDouble>(pressure);
7425 xDelete<IssmDouble>(vel);
7426 xDelete<IssmDouble>(vz);
7427 xDelete<IssmDouble>(vy);
7428 xDelete<IssmDouble>(vx);
7429 xDelete<IssmDouble>(xyz_list);
7430 xDelete<IssmDouble>(SSAvalues);
7431 xDelete<IssmDouble>(HOvalues);
7432 xDelete<int>(SSAdoflist);
7433 xDelete<int>(HOdoflist);
The documentation for this class was generated from the following files:
void GetLFSSSA(IssmDouble *L, Element *element, Gauss *gauss)
@ BasalforcingsPicoBasinIdEnum
@ FrictionCoefficientEnum
ElementVector * CreatePVectorCouplingHOFSViscous(Element *element)
@ BasalforcingsOceanTempEnum
void TransformStiffnessMatrixCoord(ElementMatrix *Ke, int cs_enum)
@ StressbalanceShelfDampeningEnum
void InputUpdateFromSolutionHO(IssmDouble *solution, Element *element)
virtual int GetElementType(void)=0
void StaticCondensation(ElementMatrix *Ke, int numindices, int *indices)
@ StressbalancePenaltyFactorEnum
ElementVector * CreatePVectorFSViscousLA(Element *element)
void GetInputListOnNodes(IssmDouble *pvalue, int enumtype)
ElementVector * CreatePVectorSSA(Element *element)
ElementVector * CreatePVectorHODrivingStress(Element *element)
ElementMatrix * CreateKMatrixFSViscous(Element *element)
ElementMatrix * CreateKMatrixSSA(Element *element)
void TransformSolutionCoord(IssmDouble *solution, int cs_enum)
virtual void NormalSection(IssmDouble *normal, IssmDouble *xyz_list)=0
void extrudefrombase_core(FemModel *femmodel)
void GetDofList(int **pdoflist, int approximation_enum, int setenum)
virtual int GetNumberOfNodes(void)=0
void InputUpdateFromSolutionL1L2(IssmDouble *solution, Element *element)
void GetBFSprimeUzawa(IssmDouble *B, Element *element, int dim, IssmDouble *xyz_list, Gauss *gauss)
#define _printf0_(StreamArgs)
@ FrictionEffectivePressureEnum
void FindParam(bool *pvalue, int paramenum)
void PetscOptionsDetermineSolverType(int *psolver_type)
ElementMatrix * CreateKMatrixSSA3d(Element *element)
ElementVector * CreatePVectorSSAFront(Element *element)
virtual void NodalFunctionsTensor(IssmDouble *basis, Gauss *gauss)=0
@ StressbalanceAnalysisEnum
virtual Gauss * NewGaussBase(int order)=0
ElementMatrix * CreateKMatrixFSFriction(Element *element)
int AddObject(Object *object)
void InputUpdateFromSolutionSSAFS(IssmDouble *solution, Element *element)
ElementVector * CreatePVectorHOFS(Element *element)
@ MeshVertexonboundaryEnum
@ StressbalanceFSreconditioningEnum
ElementMatrix * CreateKMatrixCouplingSSAFSFriction(Element *element)
@ TimesteppingTimeStepEnum
virtual int PressureInterpolation()=0
ElementVector * CreatePVectorCouplingHOFS(Element *element)
void TransformInvStiffnessMatrixCoord(ElementMatrix *Ke, int cs_enum)
ElementMatrix * CreateJacobianMatrixHO(Element *element)
void solutionsequence_schurcg(FemModel *femmodel)
virtual Element * SpawnBasalElement(void)=0
int IoCodeToEnumElementEquation(int enum_in)
ElementMatrix * NewElementMatrixCoupling(int number_nodes, int approximation_enum=NoneApproximationEnum)
ElementVector * CreatePVectorL1L2DrivingStress(Element *element)
ElementVector * CreatePVectorCouplingHOFSFriction(Element *element)
ElementMatrix * CreateKMatrixCouplingSSAHO(Element *element)
ElementMatrix * CreateKMatrixHOFriction(Element *element)
void InputUpdateFromSolutionFS(IssmDouble *solution, Element *element)
void Matrix4x4Solve(IssmDouble *X, IssmDouble *A, IssmDouble *B)
virtual void ViscosityHODerivativeEpsSquare(IssmDouble *pmu_prime, IssmDouble *epsilon, Gauss *gauss)=0
ElementMatrix * CreateKMatrixL1L2Friction(Element *element)
ElementVector * CreatePVectorHOFront(Element *element)
void InputUpdateFromSolutionHOFS(IssmDouble *solution, Element *element)
void TransformLoadVectorCoord(ElementVector *pe, int cs_enum)
ElementMatrix * CreateKMatrixHO(Element *element)
void GetDofListPressure(int **pdoflist, int setenum)
int TripleMultiply(IssmDouble *a, int nrowa, int ncola, int itrna, IssmDouble *b, int nrowb, int ncolb, int itrnb, IssmDouble *c, int nrowc, int ncolc, int itrnc, IssmDouble *d, int iaddd)
virtual void ViscosityHO(IssmDouble *pviscosity, int dim, IssmDouble *xyz_list, Gauss *gauss, Input2 *vx_input, Input2 *vy_input)=0
@ MismipFloatingMeltRateEnum
void GetBSSAHO(IssmDouble *B, Element *element, IssmDouble *xyz_list, Gauss *gauss)
ElementMatrix * CreateKMatrixCouplingHOFS(Element *element)
IssmDouble GetZcoord(IssmDouble *xyz_list, Gauss *gauss)
virtual void GaussNode(int finitelement, int iv)=0
virtual void JacobianDeterminantBase(IssmDouble *Jdet, IssmDouble *xyz_list_base, Gauss *gauss)=0
virtual void NormalBase(IssmDouble *normal, IssmDouble *xyz_list)=0
ElementVector * CreatePVectorL1L2Front(Element *element)
virtual void ViscosityFSDerivativeEpsSquare(IssmDouble *pmu_prime, IssmDouble *epsilon, Gauss *gauss)=0
@ BasalforcingsUpperwaterElevationEnum
virtual void NodalFunctionsP1Derivatives(IssmDouble *dbasis, IssmDouble *xyz_list, Gauss *gauss)=0
@ StressbalanceRiftPenaltyThresholdEnum
ElementVector * CreatePVectorFSShelf(Element *element)
@ AugmentedLagrangianThetaEnum
ElementVector * CreatePVectorFSViscous(Element *element)
@ BasalforcingsDeepwaterMeltingRateEnum
virtual void NodalFunctionsPressure(IssmDouble *basis, Gauss *gauss)=0
void AddObject(Param *newparam)
void solutionsequence_nonlinear(FemModel *femmodel, bool conserve_loads)
void GetSolutionFromInputsFS(Vector< IssmDouble > *solution, Element *element)
void GetBHOFriction(IssmDouble *B, Element *element, int dim, IssmDouble *xyz_list, Gauss *gauss)
@ StressbalanceMaxiterEnum
void GetDofListLocalPressure(int **pdoflist, int setenum)
ElementMatrix * CreateKMatrixCouplingSSAFS(Element *element)
virtual Input2 * GetInput2(int inputenum)=0
void DeleteMaterials(void)
virtual void NodalFunctions(IssmDouble *basis, Gauss *gauss)=0
@ FrictionThresholdSpeedEnum
void solutionsequence_FScoupling_nonlinear(FemModel *femmodel, bool conserve_loads)
void solutionsequence_la(FemModel *femmodel)
void StaticCondensation(int numindices, int *indices)
ElementMatrix * CreateKMatrixFSViscousXTH(Element *element)
void GetBprimeSSAFSTria(IssmDouble *Bprime, Element *element, IssmDouble *xyz_list, Gauss *gauss)
virtual void NodalFunctionsP1(IssmDouble *basis, Gauss *gauss)=0
ElementVector * NewElementVector(int approximation_enum=NoneApproximationEnum)
@ StressbalanceIsnewtonEnum
void SetParam(bool boolean, int enum_type)
void DeleteData(int num,...)
ElementVector * CreatePVectorFS(Element *element)
Param * CopyConstantObject(const char *constant_name, int param_enum)
void GetBFSvel(IssmDouble *B, Element *element, int dim, IssmDouble *xyz_list, Gauss *gauss)
@ FrictionPseudoplasticityExponentEnum
void SetValues(int ssize, int *list, doubletype *values, InsMode mode)
virtual void GetGroundedPart(int *point1, IssmDouble *fraction1, IssmDouble *fraction2, bool *mainlyfloating)=0
virtual void AddInput2(int input_enum, IssmDouble *values, int interpolation_enum)
ElementMatrix * CreateKMatrixSSALateralFriction(Element *element)
void GetBFS(IssmDouble *B, Element *element, int dim, IssmDouble *xyz_list, Gauss *gauss)
virtual Gauss * NewGauss(void)=0
@ FrictionPressureAdjustedTemperatureEnum
void GetLprimeSSAFS(IssmDouble *Lprime, Element *element, IssmDouble *xyz_list, Gauss *gauss)
ElementMatrix * CreateKMatrixCouplingSSAHOViscous(Element *element)
ElementMatrix * CreateJacobianMatrixSSA(Element *element)
void GetBFSUzawa(IssmDouble *B, Element *element, int dim, IssmDouble *xyz_list, Gauss *gauss)
@ BasalforcingsOceanSalinityEnum
ElementMatrix * CreateKMatrixHOFS(Element *element)
const char * EnumToStringx(int enum_in)
@ LinearFloatingMeltRateEnum
@ StressbalanceAbstolEnum
void GetLSSAFS(IssmDouble *L, Element *element, Gauss *gauss)
@ BeckmannGoosseFloatingMeltRateEnum
@ MantlePlumeGeothermalFluxEnum
int NewtonSolveDnorm(IssmDouble *pdnorm, IssmDouble c1, IssmDouble c2, IssmDouble c3, IssmDouble n, IssmDouble dnorm)
ElementMatrix * CreateKMatrixSSAFriction(Element *element)
@ BasalforcingsFloatingiceMeltingRateEnum
@ StressbalanceReltolEnum
void FindConstant(bool *pvalue, const char *constant_name)
ElementMatrix * CreateKMatrixSSA3dViscous(Element *element)
void GetInput2Value(bool *pvalue, int enum_type)
virtual void ViscosityFS(IssmDouble *pviscosity, int dim, IssmDouble *xyz_list, Gauss *gauss, Input2 *vx_input, Input2 *vy_input, Input2 *vz_input)=0
void GetVerticesCoordinates(IssmDouble **xyz_list)
Node * GetNode(int nodeindex)
virtual void SetTemporaryElementType(int element_type_in)=0
void GetBFSprimevel(IssmDouble *B, Element *element, int dim, IssmDouble *xyz_list, Gauss *gauss)
@ BasalforcingsPicoOverturningCoeffEnum
void GetBFSprime(IssmDouble *B, Element *element, int dim, IssmDouble *xyz_list, Gauss *gauss)
ElementMatrix * CreateKMatrixFSShelf(Element *element)
@ SpatialLinearFloatingMeltRateEnum
ElementMatrix * CreateKMatrixCouplingSSAHOFriction(Element *element)
ElementVector * CreatePVectorHO(Element *element)
void GetBSSAFS(IssmDouble *B, Element *element, IssmDouble *xyz_list, Gauss *gauss)
virtual bool IsFaceOnBoundary(void)=0
virtual void GetIcefrontCoordinates(IssmDouble **pxyz_front, IssmDouble *xyz_list, int levelsetenum)=0
virtual int NumberofNodesPressure(void)=0
void FetchData(bool *pboolean, const char *data_name)
void GetBSSA(IssmDouble *B, Element *element, int dim, IssmDouble *xyz_list, Gauss *gauss)
void CreateNodes(Nodes *nodes, IoModel *iomodel, bool isamr=false)
ElementVector * CreatePVectorSSAFS(Element *element)
void GetBSSAprime(IssmDouble *B, Element *element, int dim, IssmDouble *xyz_list, Gauss *gauss)
virtual void SynchronizeGaussBase(Gauss *gauss)=0
@ MaterialsRhoSeawaterEnum
void InputUpdateFromSolutionSSAHO(IssmDouble *solution, Element *element)
@ BasalforcingsIsmip6Enum
ElementMatrix * CreateKMatrixHOViscous(Element *element)
ElementVector * CreatePVectorSSADrivingStress(Element *element)
@ NoFrictionOnPartiallyFloatingEnum
virtual int ObjectEnum()=0
virtual void NodalFunctionsDerivativesVelocity(IssmDouble *dbasis, IssmDouble *xyz_list, Gauss *gauss)=0
void GetBprimeSSAFS(IssmDouble *Bprime, Element *element, IssmDouble *xyz_list, Gauss *gauss)
void GetDofListLocal(int **pdoflist, int approximation_enum, int setenum)
ElementMatrix * CreateKMatrixFSViscousLA(Element *element)
ElementVector * CreatePVectorCouplingSSAFSViscous(Element *element)
@ FlowequationBorderFSEnum
virtual void ViscosityL1L2(IssmDouble *pviscosity, IssmDouble *xyz_list, Gauss *gauss, Input2 *vx_input, Input2 *vy_input, Input2 *surf)=0
@ DamageStressThresholdEnum
ElementMatrix * CreateKMatrixSSAViscous(Element *element)
virtual void AddBasalInput2(int input_enum, IssmDouble *values, int interpolation_enum)
void GetBSSAFSTria(IssmDouble *B, Element *element, IssmDouble *xyz_list, Gauss *gauss)
ElementMatrix * CreateJacobianMatrixFS(Element *element)
ElementMatrix * CreateKMatrixSSA3dFriction(Element *element)
virtual void GetLevelCoordinates(IssmDouble **pxyz_front, IssmDouble *xyz_list, int levelsetenum, IssmDouble level)=0
ElementMatrix * CreateKMatrixFS(Element *element)
virtual void JacobianDeterminantSurface(IssmDouble *Jdet, IssmDouble *xyz_list, Gauss *gauss)=0
void InputUpdateFromSolutionSSA(IssmDouble *solution, Element *element)
ElementMatrix * CreateKMatrixSSAHO(Element *element)
@ AndroidFrictionCoefficientEnum
@ FrictionSedimentCompressibilityCoefficientEnum
@ BasalforcingsIsmip6BasinIdEnum
virtual void NodalFunctionsVelocity(IssmDouble *basis, Gauss *gauss)=0
IssmDouble * Data(const char *data_name)
@ StressbalanceRequestedOutputsEnum
void StrainRateHO(IssmDouble *epsilon, IssmDouble *xyz_list, Gauss *gauss, Input2 *vx_input, Input2 *vy_input)
@ StressbalanceNumRequestedOutputsEnum
virtual void GetVerticesCoordinatesBase(IssmDouble **xyz_list)=0
@ BasalforcingsDeepwaterElevationEnum
#define _error_(StreamArgs)
@ HydrologyWatercolumnMaxEnum
void SetActiveNodesLSMx(FemModel *femmodel)
void GetCFS(IssmDouble *B, Element *element, int dim, IssmDouble *xyz_list, Gauss *gauss)
int numberofverticalfaces
@ StressbalanceRestolEnum
IssmDouble GetYcoord(IssmDouble *xyz_list, Gauss *gauss)
virtual int begin(void)=0
ElementMatrix * CreateKMatrixL1L2(Element *element)
ElementVector * CreatePVectorCouplingSSAFS(Element *element)
bool VerboseSolution(void)
Object * GetObjectByOffset(int offset)
@ MeshVertexonsurfaceEnum
virtual void ViscositySSA(IssmDouble *pviscosity, int dim, IssmDouble *xyz_list, Gauss *gauss, Input2 *vx_input, Input2 *vy_input)=0
IssmDouble min(IssmDouble a, IssmDouble b)
ElementVector * CreatePVectorL1L2(Element *element)
virtual IssmDouble GetGroundedPortion(IssmDouble *xyz_list)=0
virtual void NodalFunctionsP2(IssmDouble *basis, Gauss *gauss)=0
virtual void JacobianDeterminant(IssmDouble *Jdet, IssmDouble *xyz_list, Gauss *gauss)=0
@ FrictionEffectivePressureLimitEnum
virtual void GaussPoint(int ig)=0
virtual void Update(Inputs2 *inputs2, int index, IoModel *iomodel, int analysis_counter, int analysis_type, int finite_element)=0
@ FrictionTillFrictionAngleEnum
void SetCurrentConfiguration(int configuration_type)
virtual Element * GetBasalElement(void)=0
@ SubelementFriction1Enum
void solutionsequence_la_theta(FemModel *femmodel)
virtual int VelocityInterpolation()=0
void FindParam(bool *pinteger, int enum_type)
@ SubelementFriction2Enum
virtual int FiniteElement(void)=0
ElementVector * CreatePVectorCouplingSSAFSFriction(Element *element)
@ AugmentedLagrangianREnum
void GetSolutionFromInputsHoriz(Vector< IssmDouble > *solution, Element *element)
void GetAlpha2(IssmDouble *palpha2, Gauss *gauss)
void FetchDataToInput(Inputs2 *inputs2, Elements *elements, const char *vector_name, int input_enum)
int MatrixMultiply(IssmDouble *a, int nrowa, int ncola, int itrna, IssmDouble *b, int nrowb, int ncolb, int itrnb, IssmDouble *c, int iaddc)
ElementVector * CreatePVectorFSViscousXTH(Element *element)
void GetCFSprime(IssmDouble *B, Element *element, int dim, IssmDouble *xyz_list, Gauss *gauss)
@ FrictionCoefficientcoulombEnum
void GetBSSAFriction(IssmDouble *B, Element *element, int dim, IssmDouble *xyz_list, Gauss *gauss)
void IoModelToConstraintsx(Constraints *constraints, IoModel *iomodel, const char *spc_name, int analysis_type, int finite_element, int dof)
ElementVector * CreateDVectorFS(Element *element)
void StrainRateSSA(IssmDouble *epsilon, IssmDouble *xyz_list, Gauss *gauss, Input2 *vx_input, Input2 *vy_input)
@ GroundinglineFrictionInterpolationEnum
ElementVector * CreatePVectorFSFront(Element *element)
ElementMatrix * CreateKMatrixCouplingSSAFSViscous(Element *element)
void GetLprimeFSSSA(IssmDouble *Lprime, Element *element, IssmDouble *xyz_list, Gauss *gauss)
ElementMatrix * CreateKMatrixSSAFS(Element *element)
int IoCodeToEnumVertexEquation(int enum_in)
@ AugmentedLagrangianRlambdaEnum
void Matrix3x3Solve(IssmDouble *X, IssmDouble *A, IssmDouble *B)
virtual void ViscositySSADerivativeEpsSquare(IssmDouble *pmu_prime, IssmDouble *epsilon, Gauss *gauss)=0
virtual int GetNumberOfVertices(void)=0
void GetInputListOnVertices(IssmDouble *pvalue, int enumtype)
void solutionsequence_newton(FemModel *femmodel)
ElementVector * CreatePVectorSSAHO(Element *element)
virtual bool IsIcefront(void)=0
void GetDofListLocalVelocity(int **pdoflist, int setenum)
void GetDofListVelocity(int **pdoflist, int setenum)
ElementMatrix * CreateKMatrixL1L2Viscous(Element *element)
virtual void NodalFunctionsDerivatives(IssmDouble *dbasis, IssmDouble *xyz_list, Gauss *gauss)=0
virtual void NodalFunctionsMINIDerivatives(IssmDouble *dbasis, IssmDouble *xyz_list, Gauss *gauss)=0
virtual int NumberofNodesVelocity(void)=0
ElementMatrix * NewElementMatrix(int approximation_enum=NoneApproximationEnum)