Index: /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp (revision 16992) +++ /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp (revision 16993) @@ -813,5 +813,19 @@ /*Finite Element Analysis*/ ElementMatrix* StressbalanceAnalysis::CreateJacobianMatrix(Element* element){/*{{{*/ -_error_("Not implemented"); + + int approximation; + element->GetInputValue(&approximation,ApproximationEnum); + switch(approximation){ + case SSAApproximationEnum: + return CreateJacobianMatrixSSA(element); + case HOApproximationEnum: + return CreateJacobianMatrixHO(element); + case FSApproximationEnum: + return CreateJacobianMatrixFS(element); + case NoneApproximationEnum: + return NULL; + default: + _error_("Approximation "<FindParam(&meshtype,MeshTypeEnum); + switch(meshtype){ + case Mesh2DhorizontalEnum: + basalelement = element; + break; + case Mesh3DEnum: + if(!element->IsOnBed()) return NULL; + basalelement = element->SpawnBasalElement(); + break; + default: _error_("mesh "<GetNumberOfNodes(); + + /*Initialize Element matrix, vectors and Gaussian points*/ + ElementMatrix* Ke=this->CreateKMatrixSSA(element); //Initialize Jacobian with regular SSA (first part of the Gateau derivative) + IssmDouble* dbasis = xNew(2*numnodes); + + /*Retrieve all inputs and parameters*/ + element->GetVerticesCoordinates(&xyz_list); + Input* thickness_input = basalelement->GetInput(ThicknessEnum);_assert_(thickness_input); + Input* vx_input = basalelement->GetInput(VxEnum); _assert_(vx_input); + Input* vy_input = basalelement->GetInput(VyEnum); _assert_(vy_input); + + /* Start looping on the number of gaussian points: */ + Gauss* gauss = basalelement->NewGauss(2); + for(int ig=gauss->begin();igend();ig++){ + gauss->GaussPoint(ig); + + basalelement->JacobianDeterminant(&Jdet,xyz_list,gauss); + basalelement->NodalFunctionsDerivatives(dbasis,xyz_list,gauss); + + thickness_input->GetInputValue(&thickness, gauss); + basalelement->StrainRateSSA(&epsilon[0],xyz_list,gauss,vx_input,vy_input); + basalelement->ViscositySSADerivativeEpsSquare(&mu_prime,&epsilon[0]); + eps1[0]=2*epsilon[0]+epsilon[1]; eps2[0]=epsilon[2]; + eps1[1]=epsilon[2]; eps2[1]=epsilon[0]+2*epsilon[1]; + + for(int i=0;ivalues[2*numnodes*(2*i+0)+2*j+0]+=gauss->weight*Jdet*2.*mu_prime*thickness*eps1dotdphij*eps1dotdphii; + Ke->values[2*numnodes*(2*i+0)+2*j+1]+=gauss->weight*Jdet*2.*mu_prime*thickness*eps2dotdphij*eps1dotdphii; + Ke->values[2*numnodes*(2*i+1)+2*j+0]+=gauss->weight*Jdet*2.*mu_prime*thickness*eps1dotdphij*eps2dotdphii; + Ke->values[2*numnodes*(2*i+1)+2*j+1]+=gauss->weight*Jdet*2.*mu_prime*thickness*eps2dotdphij*eps2dotdphii; + } + } + } + + /*Transform Coordinate System*/ + basalelement->TransformStiffnessMatrixCoord(Ke,XYEnum); + + /*Clean up and return*/ + delete gauss; + xDelete(xyz_list); + xDelete(dbasis); + + /*clean-up and return*/ + if(meshtype!=Mesh2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;}; + return Ke; + +}/*}}}*/ ElementMatrix* StressbalanceAnalysis::CreateKMatrixSSA(Element* element){/*{{{*/ @@ -1778,4 +1875,65 @@ /*HO*/ +ElementMatrix* StressbalanceAnalysis::CreateJacobianMatrixHO(Element* element){/*{{{*/ + + /*Intermediaries */ + IssmDouble Jdet; + IssmDouble eps1dotdphii,eps1dotdphij; + IssmDouble eps2dotdphii,eps2dotdphij; + IssmDouble mu_prime; + IssmDouble epsilon[5]; /* epsilon=[exx,eyy,exy,exz,eyz];*/ + IssmDouble eps1[3],eps2[3]; + IssmDouble *xyz_list = NULL; + + /*Fetch number of nodes and dof for this finite element*/ + int numnodes = element->GetNumberOfNodes(); + + /*Initialize Element matrix, vectors and Gaussian points*/ + ElementMatrix* Ke=this->CreateKMatrixHO(element); //Initialize Jacobian with regular HO (first part of the Gateau derivative) + IssmDouble* dbasis = xNew(3*numnodes); + + /*Retrieve all inputs and parameters*/ + element->GetVerticesCoordinates(&xyz_list); + Input* vx_input = element->GetInput(VxEnum); _assert_(vx_input); + Input* vy_input = element->GetInput(VyEnum); _assert_(vy_input); + + /* Start looping on the number of gaussian points: */ + Gauss* gauss = element->NewGauss(5); + for(int ig=gauss->begin();igend();ig++){ + gauss->GaussPoint(ig); + + element->JacobianDeterminant(&Jdet,xyz_list,gauss); + element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss); + + element->StrainRateHO(&epsilon[0],xyz_list,gauss,vx_input,vy_input); + element->ViscosityHODerivativeEpsSquare(&mu_prime,&epsilon[0]); + eps1[0]=2*epsilon[0]+epsilon[1]; eps2[0]=epsilon[2]; + eps1[1]=epsilon[2]; eps2[1]=epsilon[0]+2*epsilon[1]; + eps1[2]=epsilon[3]; eps2[2]=epsilon[4]; + + for(int i=0;ivalues[2*numnodes*(2*i+0)+2*j+0]+=gauss->weight*Jdet*2.*mu_prime*eps1dotdphij*eps1dotdphii; + Ke->values[2*numnodes*(2*i+0)+2*j+1]+=gauss->weight*Jdet*2.*mu_prime*eps2dotdphij*eps1dotdphii; + Ke->values[2*numnodes*(2*i+1)+2*j+0]+=gauss->weight*Jdet*2.*mu_prime*eps1dotdphij*eps2dotdphii; + Ke->values[2*numnodes*(2*i+1)+2*j+1]+=gauss->weight*Jdet*2.*mu_prime*eps2dotdphij*eps2dotdphii; + } + } + } + + /*Transform Coordinate System*/ + element->TransformStiffnessMatrixCoord(Ke,XYEnum); + + /*Clean up and return*/ + delete gauss; + xDelete(xyz_list); + xDelete(dbasis); + return Ke; +}/*}}}*/ ElementMatrix* StressbalanceAnalysis::CreateKMatrixHO(Element* element){/*{{{*/ @@ -2219,4 +2377,89 @@ /*FS*/ +ElementMatrix* StressbalanceAnalysis::CreateJacobianMatrixFS(Element* element){/*{{{*/ + + /*Intermediaries */ + int i,j; + IssmDouble Jdet; + IssmDouble eps1dotdphii,eps1dotdphij; + IssmDouble eps2dotdphii,eps2dotdphij; + IssmDouble eps3dotdphii,eps3dotdphij; + IssmDouble mu_prime; + IssmDouble epsilon[6]; /* epsilon=[exx,eyy,ezz,exy,exz,eyz];*/ + IssmDouble eps1[3],eps2[3],eps3[3]; + IssmDouble *xyz_list = NULL; + + /*Fetch number of nodes and dof for this finite element*/ + int vnumnodes = element->GetNumberOfNodesVelocity(); + int pnumnodes = element->GetNumberOfNodesPressure(); + int numdof = vnumnodes*NDOF3 + pnumnodes*NDOF1; + + /*Prepare coordinate system list*/ + int* cs_list = xNew(vnumnodes+pnumnodes); + for(i=0;iCreateKMatrixFS(element); //Initialize Jacobian with regular FS (first part of the Gateau derivative) + IssmDouble* dbasis = xNew(3*vnumnodes); + + /*Retrieve all inputs and parameters*/ + element->GetVerticesCoordinates(&xyz_list); + Input* vx_input = element->GetInput(VxEnum); _assert_(vx_input); + Input* vy_input = element->GetInput(VyEnum); _assert_(vy_input); + Input* vz_input = element->GetInput(VzEnum); _assert_(vz_input); + + /* Start looping on the number of gaussian points: */ + Gauss* gauss = element->NewGauss(5); + for(int ig=gauss->begin();igend();ig++){ + gauss->GaussPoint(ig); + + element->JacobianDeterminant(&Jdet,xyz_list,gauss); + element->NodalFunctionsDerivativesVelocity(dbasis,xyz_list,gauss); + + //element->StrainRateFS(&epsilon[0],xyz_list,gauss,vx_input,vy_input,vz_input); + //eps1[0]=epsilon[0]; eps2[0]=epsilon[3]; eps3[0]=epsilon[4]; + //eps1[1]=epsilon[3]; eps2[1]=epsilon[1]; eps3[1]=epsilon[5]; + //eps1[2]=epsilon[4]; eps2[2]=epsilon[5]; eps3[2]=epsilon[2]; + element->StrainRateHO(&epsilon[0],xyz_list,gauss,vx_input,vy_input); + eps1[0]=epsilon[0]; eps2[0]=epsilon[2]; eps3[0]=epsilon[3]; + eps1[1]=epsilon[2]; eps2[1]=epsilon[1]; eps3[1]=epsilon[4]; + eps1[2]=epsilon[3]; eps2[2]=epsilon[4]; eps3[2]= -epsilon[0] -epsilon[1]; + element->ViscosityFSDerivativeEpsSquare(&mu_prime,&epsilon[0]); + + for(i=0;ivalues[numdof*(3*i+0)+3*j+0]+=gauss->weight*Jdet*2*mu_prime*eps1dotdphij*eps1dotdphii; + Ke->values[numdof*(3*i+0)+3*j+1]+=gauss->weight*Jdet*2*mu_prime*eps2dotdphij*eps1dotdphii; + Ke->values[numdof*(3*i+0)+3*j+2]+=gauss->weight*Jdet*2*mu_prime*eps3dotdphij*eps1dotdphii; + + Ke->values[numdof*(3*i+1)+3*j+0]+=gauss->weight*Jdet*2*mu_prime*eps1dotdphij*eps2dotdphii; + Ke->values[numdof*(3*i+1)+3*j+1]+=gauss->weight*Jdet*2*mu_prime*eps2dotdphij*eps2dotdphii; + Ke->values[numdof*(3*i+1)+3*j+2]+=gauss->weight*Jdet*2*mu_prime*eps3dotdphij*eps2dotdphii; + + Ke->values[numdof*(3*i+2)+3*j+0]+=gauss->weight*Jdet*2*mu_prime*eps1dotdphij*eps3dotdphii; + Ke->values[numdof*(3*i+2)+3*j+1]+=gauss->weight*Jdet*2*mu_prime*eps2dotdphij*eps3dotdphii; + Ke->values[numdof*(3*i+2)+3*j+2]+=gauss->weight*Jdet*2*mu_prime*eps3dotdphij*eps3dotdphii; + } + } + } + + /*Transform Coordinate System*/ + element->TransformStiffnessMatrixCoord(Ke,cs_list); + + /*Clean up and return*/ + delete gauss; + xDelete(xyz_list); + xDelete(dbasis); + xDelete(cs_list); + return Ke; +}/*}}}*/ ElementMatrix* StressbalanceAnalysis::CreateKMatrixFS(Element* element){/*{{{*/ Index: /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.h =================================================================== --- /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.h (revision 16992) +++ /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.h (revision 16993) @@ -29,4 +29,5 @@ /*SSA*/ + ElementMatrix* CreateJacobianMatrixSSA(Element* element); ElementMatrix* CreateKMatrixSSA(Element* element); ElementMatrix* CreateKMatrixSSAViscous(Element* element); @@ -48,4 +49,5 @@ void InputUpdateFromSolutionL1L2(IssmDouble* solution,Element* element); /*HO*/ + ElementMatrix* CreateJacobianMatrixHO(Element* element); ElementMatrix* CreateKMatrixHO(Element* element); ElementMatrix* CreateKMatrixHOViscous(Element* element); @@ -59,4 +61,5 @@ void InputUpdateFromSolutionHO(IssmDouble* solution,Element* element); /*FS*/ + ElementMatrix* CreateJacobianMatrixFS(Element* element); ElementMatrix* CreateKMatrixFS(Element* element); ElementMatrix* CreateKMatrixFSViscous(Element* element); Index: /issm/trunk-jpl/src/c/classes/Elements/Element.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Element.h (revision 16992) +++ /issm/trunk-jpl/src/c/classes/Elements/Element.h (revision 16993) @@ -104,7 +104,5 @@ virtual void SetCurrentConfiguration(Elements* elements,Loads* loads,Nodes* nodes,Materials* materials,Parameters* parameters)=0; virtual void SetwiseNodeConnectivity(int* d_nz,int* o_nz,Node* node,bool* flags,int* flagsindices,int set1_enum,int set2_enum)=0; - virtual ElementMatrix* CreateKMatrix(void)=0; virtual void CreateDVector(Vector* df)=0; - virtual void CreateJacobianMatrix(Matrix* Jff)=0; virtual void ElementSizes(IssmDouble* phx,IssmDouble* phy,IssmDouble* phz)=0; virtual void EnthalpyToThermal(IssmDouble* ptemperature,IssmDouble* pwaterfraction,IssmDouble enthalpy,IssmDouble pressure)=0; Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp (revision 16992) +++ /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp (revision 16993) @@ -379,141 +379,4 @@ } /*}}}*/ -/*FUNCTION Penta::CreateKMatrix(void){{{*/ -ElementMatrix* Penta::CreateKMatrix(void){ - - /*retrieve parameters: */ - int analysis_type; - parameters->FindParam(&analysis_type,AnalysisTypeEnum); - - /*Checks in debugging {{{*/ - _assert_(this->nodes && this->material && this->matpar && this->verticalneighbors && this->parameters && this->inputs); - /*}}}*/ - - /*Skip if water element*/ - if(NoIceInElement()) return NULL; - - /*Just branch to the correct element stiffness matrix generator, according to the type of analysis we are carrying out: */ - switch(analysis_type){ - #ifdef _HAVE_STRESSBALANCE_ - case StressbalanceAnalysisEnum: - return CreateKMatrixStressbalanceHoriz(); - break; - case AdjointHorizAnalysisEnum: - return CreateKMatrixAdjointHoriz(); - break; - case StressbalanceSIAAnalysisEnum: - return CreateKMatrixStressbalanceSIA(); - break; - case StressbalanceVerticalAnalysisEnum: - return CreateKMatrixStressbalanceVert(); - break; - #endif - case L2ProjectionBaseAnalysisEnum: - return CreateBasalMassMatrix(); - break; - case MasstransportAnalysisEnum: - return CreateKMatrixMasstransport(); - break; - case FreeSurfaceTopAnalysisEnum: - return CreateKMatrixFreeSurfaceTop(); - break; - case FreeSurfaceBaseAnalysisEnum: - return CreateKMatrixFreeSurfaceBase(); - break; - #ifdef _HAVE_BALANCED_ - case BalancethicknessAnalysisEnum: - return CreateKMatrixBalancethickness(); - break; - #endif - #ifdef _HAVE_THERMAL_ - case ThermalAnalysisEnum: - return CreateKMatrixThermal(); - break; - case EnthalpyAnalysisEnum: - return CreateKMatrixEnthalpy(); - break; - case MeltingAnalysisEnum: - return CreateKMatrixMelting(); - break; - #endif - #ifdef _HAVE_HYDROLOGY_ - case HydrologyDCInefficientAnalysisEnum: - return CreateKMatrixHydrologyDCInefficient(); - break; - case HydrologyDCEfficientAnalysisEnum: - return CreateKMatrixHydrologyDCEfficient(); - break; - case L2ProjectionEPLAnalysisEnum: - return CreateEPLDomainMassMatrix(); - break; - #endif - default: - _error_("analysis " << EnumToStringx(analysis_type) << " not supported yet"); - } - - /*Make compiler happy*/ - return NULL; -} -/*}}}*/ -/*FUNCTION Penta::CreateKMatrixMasstransport {{{*/ -ElementMatrix* Penta::CreateKMatrixMasstransport(void){ - - if (!IsOnBed()) return NULL; - - /*Depth Averaging Vx and Vy*/ - this->InputDepthAverageAtBase(VxEnum,VxAverageEnum); - this->InputDepthAverageAtBase(VyEnum,VyAverageEnum); - - Tria* tria=(Tria*)SpawnTria(0); //lower face is 0, upper face is 1. - ElementMatrix* Ke=tria->CreateKMatrixMasstransport(); - delete tria->material; delete tria; - - /*Delete Vx and Vy averaged*/ - this->inputs->DeleteInput(VxAverageEnum); - this->inputs->DeleteInput(VyAverageEnum); - - /*clean up and return*/ - return Ke; -} -/*}}}*/ -/*FUNCTION Penta::CreateKMatrixFreeSurfaceTop {{{*/ -ElementMatrix* Penta::CreateKMatrixFreeSurfaceTop(void){ - - if(!IsOnSurface()) return NULL; - - Tria* tria=(Tria*)SpawnTria(1); //lower face is 0, upper face is 1. - ElementMatrix* Ke=tria->CreateKMatrixFreeSurfaceTop(); - delete tria->material; delete tria; - - /*clean up and return*/ - return Ke; -} -/*}}}*/ -/*FUNCTION Penta::CreateKMatrixFreeSurfaceBase {{{*/ -ElementMatrix* Penta::CreateKMatrixFreeSurfaceBase(void){ - - if(!IsOnBed()) return NULL; - - Tria* tria=(Tria*)SpawnTria(0); //lower face is 0, upper face is 1. - ElementMatrix* Ke=tria->CreateKMatrixFreeSurfaceBase(); - delete tria->material; delete tria; - - /*clean up and return*/ - return Ke; -} -/*}}}*/ -/*FUNCTION Penta::CreateBasalMassMatrix{{{*/ -ElementMatrix* Penta::CreateBasalMassMatrix(void){ - - if (!IsOnBed()) return NULL; - - Tria* tria=(Tria*)SpawnTria(0); //lower face is 0, upper face is 1. - ElementMatrix* Ke=tria->CreateMassMatrix(); - delete tria->material; delete tria; - - /*clean up and return*/ - return Ke; -} -/*}}}*/ /*FUNCTION Penta::CreateDVector {{{*/ void Penta::CreateDVector(Vector* df){ @@ -539,52 +402,4 @@ De->InsertIntoGlobal(df); delete De; - } -} -/*}}}*/ -/*FUNCTION Penta::CreateJacobianMatrix{{{*/ -void Penta::CreateJacobianMatrix(Matrix* Jff){ - - /*retrieve parameters: */ - ElementMatrix* Ke=NULL; - int analysis_type; - parameters->FindParam(&analysis_type,AnalysisTypeEnum); - - /*Checks in debugging {{{*/ - _assert_(this->nodes && this->material && this->matpar && this->verticalneighbors && this->parameters && this->inputs); - /*}}}*/ - - /*Skip if water element*/ - if(NoIceInElement()) return; - - /*Just branch to the correct element stiffness matrix generator, according to the type of analysis we are carrying out: */ - switch(analysis_type){ -#ifdef _HAVE_STRESSBALANCE_ - case StressbalanceAnalysisEnum: - Ke=CreateJacobianStressbalanceHoriz(); - break; -#endif - default: - _error_("analysis " << analysis_type << " (" << EnumToStringx(analysis_type) << ") not supported yet"); - } - - /*Add to global matrix*/ - if(Ke){ - /*Condense if requested*/ - if(this->element_type==MINIcondensedEnum){ - int indices[3]={18,19,20}; - Ke->StaticCondensation(3,&indices[0]); - } - else if(this->element_type==P1bubblecondensedEnum){ - int size = nodes[6]->GetNumberOfDofs(NoneApproximationEnum,GsetEnum); - int offset = 0; - for(int i=0;i<6;i++) offset+=nodes[i]->GetNumberOfDofs(NoneApproximationEnum,GsetEnum); - int* indices=xNew(size); - for(int i=0;iStaticCondensation(size,indices); - xDelete(indices); - } - - Ke->AddToGlobal(Jff); - delete Ke; } } @@ -3437,455 +3252,4 @@ #ifdef _HAVE_THERMAL_ -/*FUNCTION Penta::CreateKMatrixEnthalpy {{{*/ -ElementMatrix* Penta::CreateKMatrixEnthalpy(void){ - - /*compute all stiffness matrices for this element*/ - ElementMatrix* Ke1=CreateKMatrixEnthalpyVolume(); - ElementMatrix* Ke2=CreateKMatrixEnthalpyShelf(); - ElementMatrix* Ke =new ElementMatrix(Ke1,Ke2); - - /*clean-up and return*/ - delete Ke1; - delete Ke2; - return Ke; -} -/*}}}*/ -/*FUNCTION Penta::CreateKMatrixEnthalpyVolume {{{*/ -ElementMatrix* Penta::CreateKMatrixEnthalpyVolume(void){ - - /*Constants*/ - const int numdof=NDOF1*NUMVERTICES; - - /*Intermediaries */ - int stabilization; - int i,j; - IssmDouble Jdet,u,v,w,um,vm,wm; - IssmDouble h,hx,hy,hz,vx,vy,vz,vel; - IssmDouble gravity,rho_ice,rho_water; - IssmDouble epsvel=2.220446049250313e-16; - IssmDouble heatcapacity,thermalconductivity,dt; - IssmDouble enthalpy[NUMVERTICES], pressure[NUMVERTICES]; - IssmDouble latentheat,kappa; - IssmDouble tau_parameter,diameter; - IssmDouble xyz_list[NUMVERTICES][3]; - IssmDouble B_conduct[3][numdof]; - IssmDouble B_advec[3][numdof]; - IssmDouble Bprime_advec[3][numdof]; - IssmDouble L[numdof]; - IssmDouble dbasis[3][6]; - IssmDouble D_scalar_conduct,D_scalar_advec; - IssmDouble D_scalar_trans,D_scalar_stab; - IssmDouble D[3][3]; - IssmDouble K[3][3]={0.0}; - GaussPenta *gauss=NULL; - - /*Initialize Element matrix*/ - ElementMatrix* Ke=new ElementMatrix(nodes,NUMVERTICES,this->parameters,NoneApproximationEnum); - - /*Retrieve all inputs and parameters*/ - ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); - rho_water=matpar->GetRhoWater(); - rho_ice=matpar->GetRhoIce(); - gravity=matpar->GetG(); - heatcapacity=matpar->GetHeatCapacity(); - latentheat=matpar->GetLatentHeat(); - thermalconductivity=matpar->GetThermalConductivity(); - this->parameters->FindParam(&dt,TimesteppingTimeStepEnum); - this->parameters->FindParam(&stabilization,ThermalStabilizationEnum); - Input* vx_input=inputs->GetInput(VxEnum); _assert_(vx_input); - Input* vy_input=inputs->GetInput(VyEnum); _assert_(vy_input); - Input* vz_input=inputs->GetInput(VzEnum); _assert_(vz_input); - Input* vxm_input=inputs->GetInput(VxMeshEnum); _assert_(vxm_input); - Input* vym_input=inputs->GetInput(VyMeshEnum); _assert_(vym_input); - Input* vzm_input=inputs->GetInput(VzMeshEnum); _assert_(vzm_input); - if (stabilization==2) diameter=MinEdgeLength(&xyz_list[0][0]); - - /* Start looping on the number of gaussian points: */ - gauss=new GaussPenta(2,2); - for(int ig=gauss->begin();igend();ig++){ - - gauss->GaussPoint(ig); - - GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss); - - /*Conduction: */ - GetBConduct(&B_conduct[0][0],&xyz_list[0][0],gauss); - GetInputListOnVertices(&enthalpy[0],EnthalpyPicardEnum); - GetInputListOnVertices(&pressure[0],PressureEnum); - kappa=matpar->GetEnthalpyDiffusionParameterVolume(NUMVERTICES,enthalpy,pressure); _assert_(kappa>0.); - - D_scalar_conduct=gauss->weight*Jdet*kappa/rho_ice; - if(reCast(dt)) D_scalar_conduct=D_scalar_conduct*dt; - - D[0][0]=D_scalar_conduct; D[0][1]=0; D[0][2]=0; - D[1][0]=0; D[1][1]=D_scalar_conduct; D[1][2]=0; - D[2][0]=0; D[2][1]=0; D[2][2]=D_scalar_conduct; - - TripleMultiply(&B_conduct[0][0],3,numdof,1, - &D[0][0],3,3,0, - &B_conduct[0][0],3,numdof,0, - &Ke->values[0],1); - - /*Advection: */ - GetBAdvec(&B_advec[0][0],&xyz_list[0][0],gauss); - GetBprimeAdvec(&Bprime_advec[0][0],&xyz_list[0][0],gauss); - - vx_input->GetInputValue(&u, gauss); vxm_input->GetInputValue(&um,gauss); vx=u-um; - vy_input->GetInputValue(&v, gauss); vym_input->GetInputValue(&vm,gauss); vy=v-vm; - vz_input->GetInputValue(&w, gauss); vzm_input->GetInputValue(&wm,gauss); vz=w-wm; - - D_scalar_advec=gauss->weight*Jdet; - if(reCast(dt)) D_scalar_advec=D_scalar_advec*dt; - - D[0][0]=D_scalar_advec*vx;D[0][1]=0; D[0][2]=0; - D[1][0]=0; D[1][1]=D_scalar_advec*vy;D[1][2]=0; - D[2][0]=0; D[2][1]=0; D[2][2]=D_scalar_advec*vz; - - TripleMultiply(&B_advec[0][0],3,numdof,1, - &D[0][0],3,3,0, - &Bprime_advec[0][0],3,numdof,0, - &Ke->values[0],1); - - /*Transient: */ - if(reCast(dt)){ - GetNodalFunctionsP1(&L[0], gauss); - D_scalar_trans=gauss->weight*Jdet; - - TripleMultiply(&L[0],numdof,1,0, - &D_scalar_trans,1,1,0, - &L[0],1,numdof,0, - &Ke->values[0],1); - } - - /*Artificial diffusivity*/ - if(stabilization==1){ - /*Build K: */ - ElementSizes(&hx,&hy,&hz); - vel=sqrt(vx*vx + vy*vy + vz*vz)+1.e-14; - h=sqrt( pow(hx*vx/vel,2) + pow(hy*vy/vel,2) + pow(hz*vz/vel,2)); - - K[0][0]=h/(2*vel)*vx*vx; K[0][1]=h/(2*vel)*vx*vy; K[0][2]=h/(2*vel)*vx*vz; - K[1][0]=h/(2*vel)*vy*vx; K[1][1]=h/(2*vel)*vy*vy; K[1][2]=h/(2*vel)*vy*vz; - K[2][0]=h/(2*vel)*vz*vx; K[2][1]=h/(2*vel)*vz*vy; K[2][2]=h/(2*vel)*vz*vz; - - D_scalar_stab=gauss->weight*Jdet; - if(reCast(dt)) D_scalar_stab=D_scalar_stab*dt; - for(i=0;i<3;i++) for(j=0;j<3;j++) K[i][j] = D_scalar_stab*K[i][j]; - - GetBprimeAdvec(&Bprime_advec[0][0],&xyz_list[0][0],gauss); - - TripleMultiply(&Bprime_advec[0][0],3,numdof,1, - &K[0][0],3,3,0, - &Bprime_advec[0][0],3,numdof,0, - &Ke->values[0],1); - } - else if(stabilization==2){ - GetNodalFunctionsP1Derivatives(&dbasis[0][0],&xyz_list[0][0], gauss); - tau_parameter=StabilizationParameter(u-um,v-vm,w-wm,diameter,kappa/rho_ice); - - for(i=0;ivalues[i*numdof+j]+=tau_parameter*D_scalar_advec* - ((u-um)*dbasis[0][i]+(v-vm)*dbasis[1][i]+(w-wm)*dbasis[2][i])*((u-um)*dbasis[0][j]+(v-vm)*dbasis[1][j]+(w-wm)*dbasis[2][j]); - } - } - if(reCast(dt)){ - for(i=0;ivalues[i*numdof+j]+=tau_parameter*D_scalar_trans*L[j]*((u-um)*dbasis[0][i]+(v-vm)*dbasis[1][i]+(w-wm)*dbasis[2][i]); - } - } - } - } - } - - /*Clean up and return*/ - delete gauss; - return Ke; -} -/*}}}*/ -/*FUNCTION Penta::CreateKMatrixEnthalpyShelf {{{*/ -ElementMatrix* Penta::CreateKMatrixEnthalpyShelf(void){ - - /*Constants*/ - const int numdof=NDOF1*NUMVERTICES; - - /*Intermediaries */ - int i,j; - IssmDouble mixed_layer_capacity,thermal_exchange_velocity; - IssmDouble rho_ice,rho_water,heatcapacity; - IssmDouble Jdet2d,dt; - IssmDouble xyz_list[NUMVERTICES][3]; - IssmDouble xyz_list_tria[NUMVERTICES2D][3]; - IssmDouble basis[NUMVERTICES]; - IssmDouble D_scalar; - GaussPenta *gauss=NULL; - - /*Initialize Element matrix and return if necessary*/ - if (!IsOnBed() || !IsFloating()) return NULL; - ElementMatrix* Ke=new ElementMatrix(nodes,NUMVERTICES,this->parameters,NoneApproximationEnum); - - /*Retrieve all inputs and parameters*/ - this->parameters->FindParam(&dt,TimesteppingTimeStepEnum); - mixed_layer_capacity=matpar->GetMixedLayerCapacity(); - thermal_exchange_velocity=matpar->GetThermalExchangeVelocity(); - rho_water=matpar->GetRhoWater(); - rho_ice=matpar->GetRhoIce(); - heatcapacity=matpar->GetHeatCapacity(); - ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); - for(i=0;ibegin();igend();ig++){ - - gauss->GaussPoint(ig); - - GetTriaJacobianDeterminant(&Jdet2d, &xyz_list_tria[0][0], gauss); - GetNodalFunctionsP1(&basis[0], gauss); - - D_scalar=gauss->weight*Jdet2d*rho_water*mixed_layer_capacity*thermal_exchange_velocity/(rho_ice*heatcapacity); - if(reCast(dt)) D_scalar=dt*D_scalar; - - TripleMultiply(&basis[0],numdof,1,0, - &D_scalar,1,1,0, - &basis[0],1,numdof,0, - &Ke->values[0],1); - } - - /*Clean up and return*/ - delete gauss; - return Ke; -} -/*}}}*/ -/*FUNCTION Penta::CreateKMatrixMelting {{{*/ -ElementMatrix* Penta::CreateKMatrixMelting(void){ - - if (!IsOnBed()) return NULL; - - Tria* tria=(Tria*)SpawnTria(0); //lower face is 0, upper face is 1. - ElementMatrix* Ke=tria->CreateKMatrixMelting(); - - delete tria->material; delete tria; - return Ke; -} -/*}}}*/ -/*FUNCTION Penta::CreateKMatrixThermal {{{*/ -ElementMatrix* Penta::CreateKMatrixThermal(void){ - - /*compute all stiffness matrices for this element*/ - ElementMatrix* Ke1=CreateKMatrixThermalVolume(); - ElementMatrix* Ke2=CreateKMatrixThermalShelf(); - ElementMatrix* Ke =new ElementMatrix(Ke1,Ke2); - - /*clean-up and return*/ - delete Ke1; - delete Ke2; - return Ke; -} -/*}}}*/ -/*FUNCTION Penta::CreateKMatrixThermalVolume {{{*/ -ElementMatrix* Penta::CreateKMatrixThermalVolume(void){ - - /*Constants*/ - const int numdof=NDOF1*NUMVERTICES; - - /*Intermediaries */ - int stabilization; - int i,j; - IssmDouble Jdet,u,v,w,um,vm,wm,vel; - IssmDouble h,hx,hy,hz,vx,vy,vz; - IssmDouble gravity,rho_ice,rho_water,kappa; - IssmDouble heatcapacity,thermalconductivity,dt; - IssmDouble tau_parameter,diameter; - IssmDouble xyz_list[NUMVERTICES][3]; - IssmDouble B_conduct[3][numdof]; - IssmDouble B_advec[3][numdof]; - IssmDouble Bprime_advec[3][numdof]; - IssmDouble L[numdof]; - IssmDouble dbasis[3][6]; - IssmDouble D_scalar_conduct,D_scalar_advec; - IssmDouble D_scalar_trans,D_scalar_stab; - IssmDouble D[3][3]; - IssmDouble K[3][3]={0.0}; - GaussPenta *gauss=NULL; - - /*Initialize Element matrix*/ - ElementMatrix* Ke=new ElementMatrix(nodes,NUMVERTICES,this->parameters,NoneApproximationEnum); - - /*Retrieve all inputs and parameters*/ - ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); - rho_water=matpar->GetRhoWater(); - rho_ice=matpar->GetRhoIce(); - gravity=matpar->GetG(); - heatcapacity=matpar->GetHeatCapacity(); - thermalconductivity=matpar->GetThermalConductivity(); - kappa=thermalconductivity/(rho_ice*heatcapacity); - this->parameters->FindParam(&dt,TimesteppingTimeStepEnum); - this->parameters->FindParam(&stabilization,ThermalStabilizationEnum); - Input* vx_input=inputs->GetInput(VxEnum); _assert_(vx_input); - Input* vy_input=inputs->GetInput(VyEnum); _assert_(vy_input); - Input* vz_input=inputs->GetInput(VzEnum); _assert_(vz_input); - Input* vxm_input=inputs->GetInput(VxMeshEnum); _assert_(vxm_input); - Input* vym_input=inputs->GetInput(VyMeshEnum); _assert_(vym_input); - Input* vzm_input=inputs->GetInput(VzMeshEnum); _assert_(vzm_input); - if (stabilization==2) diameter=MinEdgeLength(&xyz_list[0][0]); - - /* Start looping on the number of gaussian points: */ - gauss=new GaussPenta(2,2); - for(int ig=gauss->begin();igend();ig++){ - - gauss->GaussPoint(ig); - - GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss); - - /*Conduction: */ - - GetBConduct(&B_conduct[0][0],&xyz_list[0][0],gauss); - - D_scalar_conduct=gauss->weight*Jdet*kappa; - if(reCast(dt)) D_scalar_conduct=D_scalar_conduct*dt; - - D[0][0]=D_scalar_conduct; D[0][1]=0; D[0][2]=0; - D[1][0]=0; D[1][1]=D_scalar_conduct; D[1][2]=0; - D[2][0]=0; D[2][1]=0; D[2][2]=D_scalar_conduct; - - TripleMultiply(&B_conduct[0][0],3,numdof,1, - &D[0][0],3,3,0, - &B_conduct[0][0],3,numdof,0, - &Ke->values[0],1); - - /*Advection: */ - GetBAdvec(&B_advec[0][0],&xyz_list[0][0],gauss); - GetBprimeAdvec(&Bprime_advec[0][0],&xyz_list[0][0],gauss); - - vx_input->GetInputValue(&u, gauss); vxm_input->GetInputValue(&um,gauss); vx=u-um; - vy_input->GetInputValue(&v, gauss); vym_input->GetInputValue(&vm,gauss); vy=v-vm; - vz_input->GetInputValue(&w, gauss); vzm_input->GetInputValue(&wm,gauss); vz=w-wm; - - D_scalar_advec=gauss->weight*Jdet; - if(reCast(dt)) D_scalar_advec=D_scalar_advec*dt; - - D[0][0]=D_scalar_advec*vx; D[0][1]=0; D[0][2]=0; - D[1][0]=0; D[1][1]=D_scalar_advec*vy; D[1][2]=0; - D[2][0]=0; D[2][1]=0; D[2][2]=D_scalar_advec*vz; - - TripleMultiply(&B_advec[0][0],3,numdof,1, - &D[0][0],3,3,0, - &Bprime_advec[0][0],3,numdof,0, - &Ke->values[0],1); - - /*Transient: */ - if(reCast(dt)){ - GetNodalFunctionsP1(&L[0], gauss); - D_scalar_trans=gauss->weight*Jdet; - - TripleMultiply(&L[0],numdof,1,0, - &D_scalar_trans,1,1,0, - &L[0],1,numdof,0, - &Ke->values[0],1); - } - - /*Artifficial diffusivity*/ - if(stabilization==1){ - /*Build K: */ - ElementSizes(&hx,&hy,&hz); - vel=sqrt(vx*vx + vy*vy + vz*vz)+1.e-14; - h=sqrt( pow(hx*vx/vel,2) + pow(hy*vy/vel,2) + pow(hz*vz/vel,2)); - - K[0][0]=h/(2*vel)*fabs(vx*vx); K[0][1]=h/(2*vel)*fabs(vx*vy); K[0][2]=h/(2*vel)*fabs(vx*vz); - K[1][0]=h/(2*vel)*fabs(vy*vx); K[1][1]=h/(2*vel)*fabs(vy*vy); K[1][2]=h/(2*vel)*fabs(vy*vz); - K[2][0]=h/(2*vel)*fabs(vz*vx); K[2][1]=h/(2*vel)*fabs(vz*vy); K[2][2]=h/(2*vel)*fabs(vz*vz); - - D_scalar_stab=gauss->weight*Jdet; - if(reCast(dt)) D_scalar_stab=D_scalar_stab*dt; - for(i=0;i<3;i++) for(j=0;j<3;j++) K[i][j] = D_scalar_stab*K[i][j]; - - GetBprimeAdvec(&Bprime_advec[0][0],&xyz_list[0][0],gauss); - - TripleMultiply(&Bprime_advec[0][0],3,numdof,1, - &K[0][0],3,3,0, - &Bprime_advec[0][0],3,numdof,0, - &Ke->values[0],1); - } - else if(stabilization==2){ - GetNodalFunctionsP1Derivatives(&dbasis[0][0],&xyz_list[0][0], gauss); - tau_parameter=StabilizationParameter(u-um,v-vm,w-wm,diameter,kappa); - - for(i=0;ivalues[i*numdof+j]+=tau_parameter*D_scalar_advec* - ((u-um)*dbasis[0][i]+(v-vm)*dbasis[1][i]+(w-wm)*dbasis[2][i])*((u-um)*dbasis[0][j]+(v-vm)*dbasis[1][j]+(w-wm)*dbasis[2][j]); - } - } - if(reCast(dt)){ - for(i=0;ivalues[i*numdof+j]+=tau_parameter*D_scalar_trans*L[j]*((u-um)*dbasis[0][i]+(v-vm)*dbasis[1][i]+(w-wm)*dbasis[2][i]); - } - } - } - } - } - - /*Clean up and return*/ - delete gauss; - return Ke; -} -/*}}}*/ -/*FUNCTION Penta::CreateKMatrixThermalShelf {{{*/ -ElementMatrix* Penta::CreateKMatrixThermalShelf(void){ - - /*Constants*/ - const int numdof=NDOF1*NUMVERTICES; - - /*Intermediaries */ - int i,j; - IssmDouble mixed_layer_capacity,thermal_exchange_velocity; - IssmDouble rho_ice,rho_water,heatcapacity; - IssmDouble Jdet2d,dt; - IssmDouble xyz_list[NUMVERTICES][3]; - IssmDouble xyz_list_tria[NUMVERTICES2D][3]; - IssmDouble basis[NUMVERTICES]; - IssmDouble D_scalar; - GaussPenta *gauss=NULL; - - /*Initialize Element matrix and return if necessary*/ - if (!IsOnBed() || !IsFloating()) return NULL; - ElementMatrix* Ke=new ElementMatrix(nodes,NUMVERTICES,this->parameters,NoneApproximationEnum); - - /*Retrieve all inputs and parameters*/ - this->parameters->FindParam(&dt,TimesteppingTimeStepEnum); - mixed_layer_capacity=matpar->GetMixedLayerCapacity(); - thermal_exchange_velocity=matpar->GetThermalExchangeVelocity(); - rho_water=matpar->GetRhoWater(); - rho_ice=matpar->GetRhoIce(); - heatcapacity=matpar->GetHeatCapacity(); - ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); - for(i=0;ibegin();igend();ig++){ - - gauss->GaussPoint(ig); - - GetTriaJacobianDeterminant(&Jdet2d, &xyz_list_tria[0][0], gauss); - GetNodalFunctionsP1(&basis[0], gauss); - - D_scalar=gauss->weight*Jdet2d*rho_water*mixed_layer_capacity*thermal_exchange_velocity/(heatcapacity*rho_ice); - if(reCast(dt)) D_scalar=dt*D_scalar; - - TripleMultiply(&basis[0],numdof,1,0, - &D_scalar,1,1,0, - &basis[0],1,numdof,0, - &Ke->values[0],1); - } - - /*Clean up and return*/ - delete gauss; - return Ke; -} -/*}}}*/ /*FUNCTION Penta::UpdateBasalConstraintsEnthalpy{{{*/ void Penta::UpdateBasalConstraintsEnthalpy(void){ @@ -4282,199 +3646,4 @@ }/*}}}*/ -/*FUNCTION Penta::CreateKMatrixAdjointHoriz{{{*/ -ElementMatrix* Penta::CreateKMatrixAdjointHoriz(void){ - - int approximation; - inputs->GetInputValue(&approximation,ApproximationEnum); - - switch(approximation){ - case SSAApproximationEnum: - return CreateKMatrixAdjointSSA2d(); - case HOApproximationEnum: - return CreateKMatrixAdjointHO(); - case FSApproximationEnum: - return CreateKMatrixAdjointFS(); - case NoneApproximationEnum: - return NULL; - default: - _error_("Approximation " << EnumToStringx(approximation) << " not supported yet"); - } -} -/*}}}*/ -/*FUNCTION Penta::CreateKMatrixAdjointSSA2d{{{*/ -ElementMatrix* Penta::CreateKMatrixAdjointSSA2d(void){ - - /*Figure out if this penta is collapsed. If so, then bailout, except if it is at the - bedrock, in which case we spawn a tria element using the 3 first nodes, and use it to build - the stiffness matrix. */ - if (!IsOnBed()) return NULL; - - /*Depth average some fields*/ - switch(this->material->ObjectEnum()){ - case MaticeEnum: - this->InputDepthAverageAtBase(MaterialsRheologyBEnum,MaterialsRheologyBbarEnum,MaterialsEnum); - this->InputDepthAverageAtBase(DamageDEnum,DamageDbarEnum,MaterialsEnum); - break; - default: - _error_("material "<material->ObjectEnum())<<" not supported"); - } - - /*Call Tria function*/ - Tria* tria=(Tria*)SpawnTria(0); //lower face is 0, upper face is 1. - ElementMatrix* Ke=tria->CreateKMatrixAdjointSSA(); - delete tria->material; delete tria; - - /*Delete averaged fields*/ - this->material->inputs->DeleteInput(MaterialsRheologyBbarEnum); - this->material->inputs->DeleteInput(DamageDbarEnum); - - /*clean up and return*/ - return Ke; -} -/*}}}*/ -/*FUNCTION Penta::CreateKMatrixAdjointHO{{{*/ -ElementMatrix* Penta::CreateKMatrixAdjointHO(void){ - - /*Intermediaries */ - int i,j; - bool incomplete_adjoint; - IssmDouble xyz_list[NUMVERTICES][3]; - IssmDouble Jdet; - IssmDouble eps1dotdphii,eps1dotdphij; - IssmDouble eps2dotdphii,eps2dotdphij; - IssmDouble mu_prime; - IssmDouble epsilon[5]; /* epsilon=[exx,eyy,exy,exz,eyz];*/ - IssmDouble eps1[3],eps2[3]; - IssmDouble dphi[3][NUMVERTICES]; - GaussPenta *gauss=NULL; - - /*Initialize Jacobian with regular HO (first part of the Gateau derivative)*/ - parameters->FindParam(&incomplete_adjoint,InversionIncompleteAdjointEnum); - ElementMatrix* Ke=CreateKMatrixStressbalanceHO(); - if(incomplete_adjoint) return Ke; - - /*Retrieve all inputs and parameters*/ - ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); - Input* vx_input=inputs->GetInput(VxEnum); _assert_(vx_input); - Input* vy_input=inputs->GetInput(VyEnum); _assert_(vy_input); - - /* Start looping on the number of gaussian points: */ - gauss=new GaussPenta(5,5); - for(int ig=gauss->begin();igend();ig++){ - - gauss->GaussPoint(ig); - - GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss); - GetNodalFunctionsP1Derivatives(&dphi[0][0],&xyz_list[0][0],gauss); - - this->StrainRateHO(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input); - material->GetViscosityDerivativeEpsSquare(&mu_prime,&epsilon[0]); - eps1[0]=2*epsilon[0]+epsilon[1]; eps2[0]=epsilon[2]; - eps1[1]=epsilon[2]; eps2[1]=epsilon[0]+2*epsilon[1]; - eps1[2]=epsilon[3]; eps2[2]=epsilon[4]; - - for(i=0;i<6;i++){ - for(j=0;j<6;j++){ - eps1dotdphii=eps1[0]*dphi[0][i]+eps1[1]*dphi[1][i]+eps1[2]*dphi[2][i]; - eps1dotdphij=eps1[0]*dphi[0][j]+eps1[1]*dphi[1][j]+eps1[2]*dphi[2][j]; - eps2dotdphii=eps2[0]*dphi[0][i]+eps2[1]*dphi[1][i]+eps2[2]*dphi[2][i]; - eps2dotdphij=eps2[0]*dphi[0][j]+eps2[1]*dphi[1][j]+eps2[2]*dphi[2][j]; - - Ke->values[12*(2*i+0)+2*j+0]+=gauss->weight*Jdet*2*mu_prime*eps1dotdphij*eps1dotdphii; - Ke->values[12*(2*i+0)+2*j+1]+=gauss->weight*Jdet*2*mu_prime*eps2dotdphij*eps1dotdphii; - Ke->values[12*(2*i+1)+2*j+0]+=gauss->weight*Jdet*2*mu_prime*eps1dotdphij*eps2dotdphii; - Ke->values[12*(2*i+1)+2*j+1]+=gauss->weight*Jdet*2*mu_prime*eps2dotdphij*eps2dotdphii; - } - } - } - - /*Transform Coordinate System*/ - ::TransformStiffnessMatrixCoord(Ke,nodes,NUMVERTICES,XYEnum); - - /*Clean up and return*/ - delete gauss; - return Ke; -} -/*}}}*/ -/*FUNCTION Penta::CreateKMatrixAdjointFS{{{*/ -ElementMatrix* Penta::CreateKMatrixAdjointFS(void){ - - /*Intermediaries */ - int i,j; - bool incomplete_adjoint; - IssmDouble xyz_list[NUMVERTICES][3]; - IssmDouble Jdet; - IssmDouble eps1dotdphii,eps1dotdphij; - IssmDouble eps2dotdphii,eps2dotdphij; - IssmDouble eps3dotdphii,eps3dotdphij; - IssmDouble mu_prime; - IssmDouble epsilon[5]; /* epsilon=[exx,eyy,exy,exz,eyz];*/ - IssmDouble eps1[3],eps2[3],eps3[3]; - IssmDouble dphi[3][NUMVERTICES]; - GaussPenta *gauss=NULL; - - /*Initialize Jacobian with regular FS (first part of the Gateau derivative)*/ - parameters->FindParam(&incomplete_adjoint,InversionIncompleteAdjointEnum); - ElementMatrix* Ke=CreateKMatrixStressbalanceFS(); - if(incomplete_adjoint) return Ke; - - /*Fetch number of nodes and dof for this finite element*/ - int vnumnodes = this->NumberofNodesVelocity(); - int pnumnodes = this->NumberofNodesPressure(); - int numdof = vnumnodes*NDOF3 + pnumnodes*NDOF1; - - /*Retrieve all inputs and parameters*/ - ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); - Input* vx_input=inputs->GetInput(VxEnum); _assert_(vx_input); - Input* vy_input=inputs->GetInput(VyEnum); _assert_(vy_input); - Input* vz_input=inputs->GetInput(VzEnum); _assert_(vz_input); - - /* Start looping on the number of gaussian points: */ - gauss=new GaussPenta(5,5); - for(int ig=gauss->begin();igend();ig++){ - - gauss->GaussPoint(ig); - - GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss); - GetNodalFunctionsP1Derivatives(&dphi[0][0],&xyz_list[0][0],gauss); - - this->StrainRateHO(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input); - material->GetViscosityDerivativeEpsSquare(&mu_prime,&epsilon[0]); - eps1[0]=epsilon[0]; eps2[0]=epsilon[2]; eps3[0]=epsilon[3]; - eps1[1]=epsilon[2]; eps2[1]=epsilon[1]; eps3[1]=epsilon[4]; - eps1[2]=epsilon[3]; eps2[2]=epsilon[4]; eps3[2]= -epsilon[0] -epsilon[1]; - - for(i=0;i<6;i++){ - for(j=0;j<6;j++){ - eps1dotdphii=eps1[0]*dphi[0][i]+eps1[1]*dphi[1][i]+eps1[2]*dphi[2][i]; - eps1dotdphij=eps1[0]*dphi[0][j]+eps1[1]*dphi[1][j]+eps1[2]*dphi[2][j]; - eps2dotdphii=eps2[0]*dphi[0][i]+eps2[1]*dphi[1][i]+eps2[2]*dphi[2][i]; - eps2dotdphij=eps2[0]*dphi[0][j]+eps2[1]*dphi[1][j]+eps2[2]*dphi[2][j]; - eps3dotdphii=eps3[0]*dphi[0][i]+eps3[1]*dphi[1][i]+eps3[2]*dphi[2][i]; - eps3dotdphij=eps3[0]*dphi[0][j]+eps3[1]*dphi[1][j]+eps3[2]*dphi[2][j]; - - Ke->values[numdof*(4*i+0)+4*j+0]+=gauss->weight*Jdet*2*mu_prime*eps1dotdphij*eps1dotdphii; - Ke->values[numdof*(4*i+0)+4*j+1]+=gauss->weight*Jdet*2*mu_prime*eps2dotdphij*eps1dotdphii; - Ke->values[numdof*(4*i+0)+4*j+2]+=gauss->weight*Jdet*2*mu_prime*eps3dotdphij*eps1dotdphii; - - Ke->values[numdof*(4*i+1)+4*j+0]+=gauss->weight*Jdet*2*mu_prime*eps1dotdphij*eps2dotdphii; - Ke->values[numdof*(4*i+1)+4*j+1]+=gauss->weight*Jdet*2*mu_prime*eps2dotdphij*eps2dotdphii; - Ke->values[numdof*(4*i+1)+4*j+2]+=gauss->weight*Jdet*2*mu_prime*eps3dotdphij*eps2dotdphii; - - Ke->values[numdof*(4*i+2)+4*j+0]+=gauss->weight*Jdet*2*mu_prime*eps1dotdphij*eps3dotdphii; - Ke->values[numdof*(4*i+2)+4*j+1]+=gauss->weight*Jdet*2*mu_prime*eps2dotdphij*eps3dotdphii; - Ke->values[numdof*(4*i+2)+4*j+2]+=gauss->weight*Jdet*2*mu_prime*eps3dotdphij*eps3dotdphii; - } - } - } - - /*Transform Coordinate System*/ - ::TransformStiffnessMatrixCoord(Ke,nodes,NUMVERTICES,XYZEnum); - - /*Clean up and return*/ - delete gauss; - return Ke; -} -/*}}}*/ /*FUNCTION Penta::GradientIndexing{{{*/ void Penta::GradientIndexing(int* indexing,int control_index){ @@ -5351,1104 +4520,4 @@ } /*}}}*/ -/*FUNCTION Penta::CreateKMatrixCouplingSSAHO{{{*/ -ElementMatrix* Penta::CreateKMatrixCouplingSSAHO(void){ - - /*compute all stiffness matrices for this element*/ - ElementMatrix* Ke1=CreateKMatrixCouplingSSAHOViscous(); - ElementMatrix* Ke2=CreateKMatrixCouplingSSAHOFriction(); - ElementMatrix* Ke=new ElementMatrix(Ke1,Ke2); - - /*clean-up and return*/ - delete Ke1; - delete Ke2; - return Ke; -} -/*}}}*/ -/*FUNCTION Penta::CreateKMatrixCouplingSSAHOViscous{{{*/ -ElementMatrix* Penta::CreateKMatrixCouplingSSAHOViscous(void){ - - /*Constants*/ - const int numnodes = 2 *NUMVERTICES; - const int numdofm = NDOF2 *NUMVERTICES2D; - const int numdofp = NDOF2 *NUMVERTICES; - const int numdoftotal = 2 *NDOF2*NUMVERTICES; - - /*Intermediaries */ - int i,j; - IssmDouble Jdet; - IssmDouble viscosity,oldviscosity,newviscosity,viscosity_overshoot; //viscosity - IssmDouble epsilon[5],oldepsilon[5]; /* epsilon=[exx,eyy,exy,exz,eyz];*/ - IssmDouble xyz_list[NUMVERTICES][3]; - IssmDouble B[3][numdofp]; - IssmDouble Bprime[3][numdofm]; - IssmDouble D[3][3]={0.0}; // material matrix, simple scalar matrix. - IssmDouble D_scalar; - IssmDouble Ke_gg[numdofp][numdofm]={0.0}; //local element stiffness matrix - IssmDouble Ke_gg_gaussian[numdofp][numdofm]; //stiffness matrix evaluated at the gaussian point. - GaussPenta *gauss=NULL; - GaussTria *gauss_tria=NULL; - Node *node_list[numnodes]; - int cs_list[numnodes]; - - /*Find penta on bed as HO must be coupled to the dofs on the bed: */ - Penta* pentabase=(Penta*)GetBasalElement(); - Tria* tria=pentabase->SpawnTria(0); //lower face is 0, upper face is 1. - - /*prepare node list*/ - for(i=0;inodes[i]; - node_list[i+1*NUMVERTICES] = this->nodes[i]; - cs_list[i+0*NUMVERTICES] = XYEnum; - cs_list[i+1*NUMVERTICES] = XYEnum; - } - - /*Initialize Element matrix*/ - ElementMatrix* Ke1=new ElementMatrix(pentabase->nodes,NUMVERTICES,this->parameters,SSAApproximationEnum); - ElementMatrix* Ke2=new ElementMatrix(this->nodes ,NUMVERTICES,this->parameters,HOApproximationEnum); - ElementMatrix* Ke =new ElementMatrix(Ke1,Ke2); - delete Ke1; delete Ke2; - - /* Get node coordinates and dof list: */ - ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); - this->parameters->FindParam(&viscosity_overshoot,StressbalanceViscosityOvershootEnum); - Input* vx_input=inputs->GetInput(VxEnum); _assert_(vx_input); - Input* vy_input=inputs->GetInput(VyEnum); _assert_(vy_input); - Input* vxold_input=inputs->GetInput(VxPicardEnum); _assert_(vxold_input); - Input* vyold_input=inputs->GetInput(VyPicardEnum); _assert_(vyold_input); - - /* Start looping on the number of gaussian points: */ - gauss=new GaussPenta(5,5); - gauss_tria=new GaussTria(); - for(int ig=gauss->begin();igend();ig++){ - - gauss->GaussPoint(ig); - gauss->SynchronizeGaussBase(gauss_tria); - - GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss); - GetBSSAHO(&B[0][0], &xyz_list[0][0], gauss); - tria->GetBprimeSSA(&Bprime[0][0], &xyz_list[0][0], gauss_tria); - - this->StrainRateHO(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input); - this->StrainRateHO(&oldepsilon[0],&xyz_list[0][0],gauss,vxold_input,vyold_input); - material->GetViscosity3d(&viscosity, &epsilon[0]); - material->GetViscosity3d(&oldviscosity, &oldepsilon[0]); - - newviscosity=viscosity+viscosity_overshoot*(viscosity-oldviscosity); - D_scalar=2*newviscosity*gauss->weight*Jdet; - for (i=0;i<3;i++) D[i][i]=D_scalar; - - TripleMultiply( &B[0][0],3,numdofp,1, - &D[0][0],3,3,0, - &Bprime[0][0],3,numdofm,0, - &Ke_gg_gaussian[0][0],0); - - for( i=0; ivalues[(i+2*numdofm)*numdoftotal+j]+=Ke_gg[i][j]; - for(i=0;ivalues[i*numdoftotal+(j+2*numdofm)]+=Ke_gg[j][i]; - - /*Transform Coordinate System*/ - ::TransformStiffnessMatrixCoord(Ke,node_list,numnodes,cs_list); - - /*Clean-up and return*/ - delete tria->material; delete tria; - delete gauss; - delete gauss_tria; - return Ke; -} -/*}}}*/ -/*FUNCTION Penta::CreateKMatrixCouplingSSAHOFriction{{{*/ -ElementMatrix* Penta::CreateKMatrixCouplingSSAHOFriction(void){ - - /*Constants*/ - const int numnodes = 2 *NUMVERTICES; - const int numdof = NDOF2 *NUMVERTICES; - const int numdoftotal = NDOF4 *NUMVERTICES; - - /*Intermediaries */ - int i,j,analysis_type; - IssmDouble Jdet2d,alpha2; - IssmDouble xyz_list[NUMVERTICES][3]; - IssmDouble xyz_list_tria[NUMVERTICES2D][3]={0.0}; - IssmDouble L[2][numdof]; - IssmDouble DL[2][2] ={{ 0,0 },{0,0}}; //for basal drag - IssmDouble DL_scalar; - IssmDouble Ke_gg[numdof][numdof] ={0.0}; - IssmDouble Ke_gg_gaussian[numdof][numdof]; //stiffness matrix contribution from drag - Friction *friction = NULL; - GaussPenta *gauss=NULL; - Node *node_list[numnodes]; - int cs_list[numnodes]; - - /*Initialize Element matrix and return if necessary*/ - if(IsFloating() || !IsOnBed()) return NULL; - ElementMatrix* Ke1=new ElementMatrix(nodes,NUMVERTICES,this->parameters,SSAApproximationEnum); - ElementMatrix* Ke2=new ElementMatrix(nodes,NUMVERTICES,this->parameters,HOApproximationEnum); - ElementMatrix* Ke=new ElementMatrix(Ke1,Ke2); - delete Ke1; delete Ke2; - - /*Prepare node list*/ - for(i=0;inodes[i]; - node_list[i+1*NUMVERTICES] = this->nodes[i]; - cs_list[i+0*NUMVERTICES] = XYEnum; - cs_list[i+1*NUMVERTICES] = XYEnum; - } - - /*retrieve inputs :*/ - ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); - for(i=0;iFindParam(&analysis_type,AnalysisTypeEnum); - Input* vx_input=inputs->GetInput(VxEnum); _assert_(vx_input); - Input* vy_input=inputs->GetInput(VyEnum); _assert_(vy_input); - Input* vz_input=inputs->GetInput(VzEnum); _assert_(vz_input); - - /*build friction object, used later on: */ - friction=new Friction(this,2); - - /* Start looping on the number of gaussian points: */ - gauss=new GaussPenta(0,1,2,2); - for(int ig=gauss->begin();igend();ig++){ - - gauss->GaussPoint(ig); - - /*Friction: */ - friction->GetAlpha2(&alpha2,gauss,vx_input,vy_input,vz_input); - - GetTriaJacobianDeterminant(&Jdet2d, &xyz_list_tria[0][0],gauss); - GetBHOFriction(&L[0][0],gauss); - - DL_scalar=alpha2*gauss->weight*Jdet2d; - for (i=0;i<2;i++) DL[i][i]=DL_scalar; - - /* Do the triple producte tL*D*L: */ - TripleMultiply( &L[0][0],2,numdof,1, - &DL[0][0],2,2,0, - &L[0][0],2,numdof,0, - &Ke_gg_gaussian[0][0],0); - - for(i=0;ivalues[i*numdoftotal+(numdof+j)]+=Ke_gg[i][j]; - for(i=0;ivalues[(i+numdof)*numdoftotal+j]+=Ke_gg[i][j]; - - /*Transform Coordinate System*/ - ::TransformStiffnessMatrixCoord(Ke,node_list,numnodes,cs_list); - - /*Clean up and return*/ - delete gauss; - delete friction; - return Ke; -} -/*}}}*/ -/*FUNCTION Penta::CreateKMatrixCouplingSSAFS{{{*/ -ElementMatrix* Penta::CreateKMatrixCouplingSSAFS(void){ - - /*compute all stiffness matrices for this element*/ - ElementMatrix* Ke1=CreateKMatrixCouplingSSAFSViscous(); - ElementMatrix* Ke2=CreateKMatrixCouplingSSAFSFriction(); - ElementMatrix* Ke =new ElementMatrix(Ke1,Ke2); - - /*clean-up and return*/ - delete Ke1; - delete Ke2; - return Ke; -} -/*}}}*/ -/*FUNCTION Penta::CreateKMatrixCouplingSSAFSViscous{{{*/ -ElementMatrix* Penta::CreateKMatrixCouplingSSAFSViscous(void){ - - /*Constants*/ - const int numdofm = NDOF2 *NUMVERTICES2D; - const int numdofs = NDOF4 *NUMVERTICES + NDOF3; - const int numdoftotal = 2 *numdofm+numdofs; - - /*Intermediaries */ - int i,j; - IssmDouble Jdet; - IssmDouble viscosity,FSreconditioning; //viscosity - IssmDouble epsilon[6]; /* epsilon=[exx,eyy,exy,exz,eyz];*/ - IssmDouble xyz_list[NUMVERTICES][3]; - IssmDouble B[4][numdofs]; - IssmDouble Bprime[4][numdofm]; - IssmDouble B2[3][numdofm]; - IssmDouble Bprime2[3][numdofs]; - IssmDouble D[4][4]={0.0}; // material matrix, simple scalar matrix. - IssmDouble D2[3][3]={0.0}; // material matrix, simple scalar matrix. - IssmDouble D_scalar; - IssmDouble Ke_gg[numdofs][numdofm]={0.0}; //local element stiffness matrix - IssmDouble Ke_gg2[numdofm][numdofs]={0.0}; //local element stiffness matrix - IssmDouble Ke_gg_gaussian[numdofs][numdofm]; //stiffness matrix evaluated at the gaussian point. - IssmDouble Ke_gg_gaussian2[numdofm][numdofs]; //stiffness matrix evaluated at the gaussian point. - GaussPenta *gauss=NULL; - GaussTria *gauss_tria=NULL; - Node *node_list[20]; - - /*Find penta on bed as FS must be coupled to the dofs on the bed: */ - Penta* pentabase=(Penta*)GetBasalElement(); - Tria* tria=pentabase->SpawnTria(0); //lower face is 0, upper face is 1. - - int vnumnodes = this->NumberofNodesVelocity(); - int pnumnodes = this->NumberofNodesPressure(); - int numnodes = 2*vnumnodes-1+pnumnodes; - - /*Prepare node list*/ - int* cs_list = xNew(2*vnumnodes+pnumnodes); - for(i=0;inodes[i]; - cs_list[i] = XYEnum; - } - for(i=0;inodes[i]; - cs_list[i+vnumnodes-1] = XYZEnum; - } - for(i=0;inodes[vnumnodes+i]; - cs_list[2*vnumnodes-1+i] = PressureEnum; - } - - /*Initialize Element matrix and return if necessary*/ - ElementMatrix* Ke1=new ElementMatrix(pentabase->nodes,NUMVERTICES, this->parameters,SSAApproximationEnum); - ElementMatrix* Ke2=new ElementMatrix(this->nodes ,2*NUMVERTICES+1,this->parameters,FSvelocityEnum); - ElementMatrix* Ke =new ElementMatrix(Ke1,Ke2); - delete Ke1; delete Ke2; - - /* Get node coordinates and dof list: */ - ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); - parameters->FindParam(&FSreconditioning,StressbalanceFSreconditioningEnum); - Input* vx_input=inputs->GetInput(VxEnum); _assert_(vx_input); - Input* vy_input=inputs->GetInput(VyEnum); _assert_(vy_input); - Input* vz_input=inputs->GetInput(VzEnum); _assert_(vz_input); - - /* Start looping on the number of gaussian points: */ - gauss=new GaussPenta(5,5); - gauss_tria=new GaussTria(); - for(int ig=gauss->begin();igend();ig++){ - - gauss->GaussPoint(ig); - gauss->SynchronizeGaussBase(gauss_tria); - - GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss); - GetBSSAFS(&B[0][0], &xyz_list[0][0], gauss); - tria->GetBprimeSSAFS(&Bprime[0][0], &xyz_list[0][0], gauss_tria); - tria->GetBSSAFS(&B2[0][0], &xyz_list[0][0], gauss_tria); - GetBprimeSSAFS(&Bprime2[0][0], &xyz_list[0][0], gauss); - - this->StrainRateFS(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input); - material->GetViscosity3dFS(&viscosity, &epsilon[0]); - - D_scalar=2*viscosity*gauss->weight*Jdet; - for (i=0;i<3;i++) D[i][i]=D_scalar; - D[3][3]=-gauss->weight*Jdet*FSreconditioning; - for (i=0;i<3;i++) D2[i][i]=D_scalar; - - TripleMultiply( &B[0][0],4,numdofs,1, - &D[0][0],4,4,0, - &Bprime[0][0],4,numdofm,0, - &Ke_gg_gaussian[0][0],0); - - TripleMultiply( &B2[0][0],3,numdofm,1, - &D2[0][0],3,3,0, - &Bprime2[0][0],3,numdofs,0, - &Ke_gg_gaussian2[0][0],0); - - for( i=0; ivalues[(i+2*numdofm)*numdoftotal+j]+=Ke_gg[i][j]; - for(i=0;ivalues[i*numdoftotal+(j+2*numdofm)]+=Ke_gg2[i][j]; - - /*Transform Coordinate System*/ - ::TransformStiffnessMatrixCoord(Ke,node_list,numnodes,cs_list); - - /*Clean-up and return*/ - xDelete(cs_list); - delete tria->material; delete tria; - delete gauss; - delete gauss_tria; - return Ke; -} -/*}}}*/ -/*FUNCTION Penta::CreateKMatrixCouplingSSAFSFriction {{{*/ -ElementMatrix* Penta::CreateKMatrixCouplingSSAFSFriction(void){ - /*Constants*/ - const int numdofs = (NUMVERTICES+1)*NDOF3 + NUMVERTICES*NDOF1; - const int numdofm = NUMVERTICES *NDOF2; - const int numdof2d = NUMVERTICES2D *NDOF3; - const int numdof2dm = NUMVERTICES2D *NDOF2; - const int numdoftot = NUMVERTICES*2 + (NUMVERTICES+1)*3 +NUMVERTICES; // HO + FS vel + FS Pressure - - /*Intermediaries */ - int i,j; - int analysis_type,approximation; - IssmDouble FSreconditioning; - IssmDouble viscosity,alpha2_gauss,Jdet2d; - IssmDouble bed_normal[3]; - IssmDouble epsilon[6]; /* epsilon=[exx,eyy,ezz,exy,exz,eyz];*/ - IssmDouble xyz_list[NUMVERTICES][3]; - IssmDouble xyz_list_tria[NUMVERTICES2D][3]; - IssmDouble LSSAFS[8][numdof2dm]; - IssmDouble LprimeSSAFS[8][numdofs]; - IssmDouble DLSSAFS[8][8]={0.0}; - IssmDouble LFSSSA[4][numdof2d]; - IssmDouble LprimeFSSSA[4][numdof2dm]; - IssmDouble DLFSSSA[4][4]={0.0}; - IssmDouble Ke_drag_gaussian[numdof2dm][numdofs]; - IssmDouble Ke_drag_gaussian2[numdof2d][numdof2dm]; - Friction* friction=NULL; - GaussPenta *gauss=NULL; - Node *node_list[20]; - - /*If on water or not FS, skip stiffness: */ - inputs->GetInputValue(&approximation,ApproximationEnum); - if(IsFloating() || !IsOnBed()) return NULL; - - int vnumnodes = this->NumberofNodesVelocity(); - int pnumnodes = this->NumberofNodesPressure(); - int numnodes = 2*vnumnodes-1+pnumnodes; - - /*Prepare node list*/ - int* cs_list = xNew(2*vnumnodes+pnumnodes); - for(i=0;inodes[i]; - cs_list[i] = XYEnum; - } - for(i=0;inodes[i]; - cs_list[i+vnumnodes-1] = XYZEnum; - } - for(i=0;inodes[vnumnodes+i]; - cs_list[2*vnumnodes-1+i] = PressureEnum; - } - - ElementMatrix* Ke1=new ElementMatrix(this->nodes,NUMVERTICES, this->parameters,SSAApproximationEnum); - ElementMatrix* Ke2=new ElementMatrix(this->nodes,vnumnodes+pnumnodes,this->parameters,FSvelocityEnum); - ElementMatrix* Ke=new ElementMatrix(Ke1,Ke2); - delete Ke1; delete Ke2; - - /*Retrieve all inputs and parameters*/ - ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); - parameters->FindParam(&analysis_type,AnalysisTypeEnum); - parameters->FindParam(&FSreconditioning,StressbalanceFSreconditioningEnum); - Input* vx_input=inputs->GetInput(VxEnum); _assert_(vx_input); - Input* vy_input=inputs->GetInput(VyEnum); _assert_(vy_input); - Input* vz_input=inputs->GetInput(VzEnum); _assert_(vz_input); - for(i=0;ibegin();igend();ig++){ - - gauss->GaussPoint(ig); - - GetTriaJacobianDeterminant(&Jdet2d, &xyz_list_tria[0][0],gauss); - GetLSSAFS(&LSSAFS[0][0], gauss); - GetLprimeSSAFS(&LprimeSSAFS[0][0], &xyz_list[0][0], gauss); - GetLFSSSA(&LFSSSA[0][0], gauss); - GetLprimeFSSSA(&LprimeFSSSA[0][0], &xyz_list[0][0], gauss); - - this->StrainRateFS(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input); - material->GetViscosity3dFS(&viscosity,&epsilon[0]); - - NormalBase(&bed_normal[0],&xyz_list_tria[0][0]); - friction->GetAlpha2(&alpha2_gauss, gauss,vx_input,vy_input,vz_input); - - DLSSAFS[0][0]=alpha2_gauss*gauss->weight*Jdet2d; - DLSSAFS[1][1]=alpha2_gauss*gauss->weight*Jdet2d; - DLSSAFS[2][2]=-alpha2_gauss*gauss->weight*Jdet2d*bed_normal[0]*bed_normal[2]; - DLSSAFS[3][3]=-alpha2_gauss*gauss->weight*Jdet2d*bed_normal[1]*bed_normal[2]; - DLSSAFS[4][4]=-2*viscosity*gauss->weight*Jdet2d*bed_normal[0]; - DLSSAFS[5][5]=-2*viscosity*gauss->weight*Jdet2d*bed_normal[1]; - DLSSAFS[6][6]=FSreconditioning*gauss->weight*Jdet2d*bed_normal[0]; - DLSSAFS[7][7]=FSreconditioning*gauss->weight*Jdet2d*bed_normal[1]; - - DLFSSSA[0][0]=alpha2_gauss*gauss->weight*Jdet2d; - DLFSSSA[1][1]=alpha2_gauss*gauss->weight*Jdet2d; - DLFSSSA[2][2]=-alpha2_gauss*gauss->weight*Jdet2d*bed_normal[0]*bed_normal[2]; - DLFSSSA[3][3]=-alpha2_gauss*gauss->weight*Jdet2d*bed_normal[1]*bed_normal[2]; - - TripleMultiply( &LSSAFS[0][0],8,numdof2dm,1, - &DLSSAFS[0][0],8,8,0, - &LprimeSSAFS[0][0],8,numdofs,0, - &Ke_drag_gaussian[0][0],0); - - TripleMultiply( &LFSSSA[0][0],4,numdof2d,1, - &DLFSSSA[0][0],4,4,0, - &LprimeFSSSA[0][0],4,numdof2dm,0, - &Ke_drag_gaussian2[0][0],0); - for(i=0;ivalues[i*numdoftot+j+numdofm]+=Ke_drag_gaussian[i][j]; - for(i=0;ivalues[(i+numdofm)*numdoftot+j]+=Ke_drag_gaussian2[i][j]; - } - - /*Transform Coordinate System*/ - ::TransformStiffnessMatrixCoord(Ke,node_list,numnodes,cs_list); - - /*Clean up and return*/ - xDelete(cs_list); - delete gauss; - delete friction; - return Ke; -} -/*}}}*/ -/*FUNCTION Penta::CreateKMatrixCouplingHOFS{{{*/ -ElementMatrix* Penta::CreateKMatrixCouplingHOFS(void){ - - /*Constants*/ - const int numnodes = 3 *NUMVERTICES+1; - const int numdofp = NDOF2 *NUMVERTICES; - const int numdofs = NDOF4 * 6 + NDOF3; - const int numdoftotal = (NDOF2+NDOF4) *NUMVERTICES + NDOF3; - - /*Intermediaries*/ - int i,j,init; - Node *node_list[NUMVERTICES*3+1]; - int cs_list[NUMVERTICES*3+1]; - int cs_list2[NUMVERTICES*2+1]; - - /*Some parameters needed*/ - init=this->element_type; - - /*prepare node list*/ - for(i=0;inodes[i]; - cs_list[i+NUMVERTICES] = XYZEnum; - cs_list2[i] = XYZEnum; - } - for(i=0;inodes[i]; - node_list[i+2*NUMVERTICES+1] = this->nodes[i+NUMVERTICES+1]; - cs_list[i] = XYEnum; - cs_list[i+2*NUMVERTICES+1] = PressureEnum; - cs_list2[i+NUMVERTICES+1] = PressureEnum; - } - - /*compute all stiffness matrices for this element*/ - ElementMatrix* Ke1=new ElementMatrix(this->nodes,NUMVERTICES,this->parameters,HOApproximationEnum); - ElementMatrix* Ke2=new ElementMatrix(this->nodes,2*NUMVERTICES+1,this->parameters,FSvelocityEnum); - ElementMatrix* Ke=new ElementMatrix(Ke1,Ke2); - delete Ke1; delete Ke2; - - /*Compute HO Matrix with P1 element type\n");*/ - this->element_type=P1Enum; - Ke1=CreateKMatrixStressbalanceHO(); TransformInvStiffnessMatrixCoord(Ke1,XYEnum); - this->element_type=init; - /*Compute FS Matrix and condense it \n");*/ - Ke2=CreateKMatrixStressbalanceFS(); TransformInvStiffnessMatrixCoord(Ke2,node_list,2*NUMVERTICES+1,cs_list2); - int indices[3]={18,19,20}; - Ke2->StaticCondensation(3,&indices[0]); - - for(i=0;ivalues[(i+numdofp)*numdoftotal+NDOF2*j+0]+=Ke2->values[i*numdofs+NDOF3*j+0]; - Ke->values[(i+numdofp)*numdoftotal+NDOF2*j+1]+=Ke2->values[i*numdofs+NDOF3*j+1]; - } - for(i=0;ivalues[i*numdoftotal+numdofp+NDOF3*j+0]+=Ke1->values[i*numdofp+NDOF2*j+0]; - Ke->values[i*numdoftotal+numdofp+NDOF3*j+1]+=Ke1->values[i*numdofp+NDOF2*j+1]; - } - - /*Transform Coordinate System*/ //Do not transform, already done in the matrices - ::TransformStiffnessMatrixCoord(Ke,node_list,numnodes,cs_list); - - /*clean-up and return*/ - delete Ke1; - delete Ke2; - return Ke; -} -//*}}}*/ -/*FUNCTION Penta::CreateKMatrixStressbalanceHoriz {{{*/ -ElementMatrix* Penta::CreateKMatrixStressbalanceHoriz(void){ - - int approximation; - inputs->GetInputValue(&approximation,ApproximationEnum); - switch(approximation){ - case SSAApproximationEnum: - return CreateKMatrixStressbalanceSSA2d(); - case L1L2ApproximationEnum: - return CreateKMatrixStressbalanceL1L2(); - case HOApproximationEnum: - return CreateKMatrixStressbalanceHO(); - case FSApproximationEnum: - return CreateKMatrixStressbalanceFS(); - case SIAApproximationEnum: - return NULL; - case NoneApproximationEnum: - return NULL; - case SSAHOApproximationEnum: - return CreateKMatrixStressbalanceSSAHO(); - case SSAFSApproximationEnum: - return CreateKMatrixStressbalanceSSAFS(); - case HOFSApproximationEnum: - return CreateKMatrixStressbalanceHOFS(); - default: - _error_("Approximation " << EnumToStringx(approximation) << " not supported yet"); - } -} -/*}}}*/ -/*FUNCTION Penta::CreateKMatrixStressbalanceSIA{{{*/ -ElementMatrix* Penta::CreateKMatrixStressbalanceSIA(void){ - - /*Intermediaries*/ - IssmDouble connectivity[2]; - IssmDouble one0,one1; - int i,i0,i1,j0,j1; - - /*Fetch number of nodes and dof for this finite element*/ - int numnodes = this->NumberofNodes(); _assert_(numnodes==6); - int numdof = numnodes*NDOF2; - - /*Initialize Element matrix*/ - ElementMatrix* Ke=new ElementMatrix(nodes,numnodes,this->parameters,NoneApproximationEnum); - - /*3 vertical edges*/ - for(i=0;i<3;i++){ - - /*2 dofs of first node*/ - i0=2*i; i1=2*i+1; - /*2 dofs of second node*/ - j0=2*(i+3); j1=2*(i+3)+1; - - /*Find connectivity for the two nodes*/ - connectivity[0]=(IssmDouble)vertices[i]->Connectivity(); - connectivity[1]=(IssmDouble)vertices[i+3]->Connectivity(); - one0=1./connectivity[0]; - one1=1./connectivity[1]; - - /*Create matrix for these two nodes*/ - if (IsOnBed() && IsOnSurface()){ - Ke->values[i0*numdof+i0] = +one0; - Ke->values[i1*numdof+i1] = +one0; - Ke->values[j0*numdof+i0] = -one1; - Ke->values[j0*numdof+j0] = +one1; - Ke->values[j1*numdof+i1] = -one1; - Ke->values[j1*numdof+j1] = +one1; - } - else if (IsOnBed()){ - Ke->values[i0*numdof+i0] = one0; - Ke->values[i1*numdof+i1] = one0; - Ke->values[j0*numdof+i0] = -2.*one1; - Ke->values[j0*numdof+j0] = +2.*one1; - Ke->values[j1*numdof+i1] = -2.*one1; - Ke->values[j1*numdof+j1] = +2.*one1; - } - else if (IsOnSurface()){ - Ke->values[j0*numdof+i0] = -one1; - Ke->values[j0*numdof+j0] = +one1; - Ke->values[j1*numdof+i1] = -one1; - Ke->values[j1*numdof+j1] = +one1; - } - else{ //node is on two horizontal layers and beams include the values only once, so the have to use half of the connectivity - Ke->values[j0*numdof+i0] = -2.*one1; - Ke->values[j0*numdof+j0] = +2.*one1; - Ke->values[j1*numdof+i1] = -2.*one1; - Ke->values[j1*numdof+j1] = +2.*one1; - } - } - - /*Clean up and return*/ - return Ke; -}/*}}}*/ -/*FUNCTION Penta::CreateKMatrixStressbalanceSSA2d{{{*/ -ElementMatrix* Penta::CreateKMatrixStressbalanceSSA2d(void){ - - /*Figure out if this penta is collapsed. If so, then bailout, except if it is at the - bedrock, in which case we spawn a tria element using the 3 first nodes, and use it to build - the stiffness matrix. */ - if (!IsOnBed()) return NULL; - - /*Depth average some fields*/ - switch(this->material->ObjectEnum()){ - case MaticeEnum: - this->InputDepthAverageAtBase(MaterialsRheologyBEnum,MaterialsRheologyBbarEnum,MaterialsEnum); - this->InputDepthAverageAtBase(DamageDEnum,DamageDbarEnum,MaterialsEnum); - break; - default: - _error_("material "<material->ObjectEnum())<<" not supported"); - } - - /*Call Tria function*/ - Tria* tria=(Tria*)SpawnTria(0); //lower face is 0, upper face is 1. - ElementMatrix* Ke=tria->CreateKMatrixStressbalanceSSA(); - delete tria->material; delete tria; - - /*Delete averaged fields*/ - this->material->inputs->DeleteInput(MaterialsRheologyBbarEnum); - this->material->inputs->DeleteInput(DamageDbarEnum); - - /*clean up and return*/ - return Ke; -} -/*}}}*/ -/*FUNCTION Penta::CreateKMatrixStressbalanceSSA3d{{{*/ -ElementMatrix* Penta::CreateKMatrixStressbalanceSSA3d(void){ - - /*compute all stiffness matrices for this element*/ - ElementMatrix* Ke1=CreateKMatrixStressbalanceSSA3dViscous(); - ElementMatrix* Ke2=CreateKMatrixStressbalanceSSA3dFriction(); - ElementMatrix* Ke =new ElementMatrix(Ke1,Ke2); - - /*clean-up and return*/ - delete Ke1; - delete Ke2; - return Ke; -} -/*}}}*/ -/*FUNCTION Penta::CreateKMatrixStressbalanceSSA3dViscous{{{*/ -ElementMatrix* Penta::CreateKMatrixStressbalanceSSA3dViscous(void){ - - /*Constants*/ - const int numdof2d=2*NUMVERTICES2D; - - /*Intermediaries */ - int i,j,approximation; - IssmDouble Jdet; - IssmDouble viscosity , oldviscosity, newviscosity, viscosity_overshoot; - IssmDouble epsilon[5],oldepsilon[5]; /* epsilon=[exx,eyy,exy,exz,eyz];*/ - IssmDouble epsilons[6]; //6 for FS - IssmDouble xyz_list[NUMVERTICES][3]; - IssmDouble B[3][numdof2d]; - IssmDouble Bprime[3][numdof2d]; - IssmDouble D[3][3]= {0.0}; // material matrix, simple scalar matrix. - IssmDouble D_scalar; - IssmDouble Ke_gg_gaussian[numdof2d][numdof2d]; - Tria *tria = NULL; - Penta *pentabase = NULL; - GaussPenta *gauss = NULL; - GaussTria *gauss_tria = NULL; - - /*Find penta on bed as this is a SSA elements: */ - pentabase=(Penta*)GetBasalElement(); - tria=pentabase->SpawnTria(0); //lower face is 0, upper face is 1. - - /*Initialize Element matrix*/ - ElementMatrix* Ke=new ElementMatrix(tria->nodes,NUMVERTICES2D,this->parameters,SSAApproximationEnum); - inputs->GetInputValue(&approximation,ApproximationEnum); - - /*Retrieve all inputs and parameters*/ - ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); - this->parameters->FindParam(&viscosity_overshoot,StressbalanceViscosityOvershootEnum); - Input* vx_input=inputs->GetInput(VxEnum); _assert_(vx_input); - Input* vy_input=inputs->GetInput(VyEnum); _assert_(vy_input); - Input* vxold_input=inputs->GetInput(VxPicardEnum); _assert_(vxold_input); - Input* vyold_input=inputs->GetInput(VyPicardEnum); _assert_(vyold_input); - Input* vz_input=inputs->GetInput(VzEnum); _assert_(vz_input); - - /* Start looping on the number of gaussian points: */ - gauss=new GaussPenta(5,5); - gauss_tria=new GaussTria(); - for(int ig=gauss->begin();igend();ig++){ - - gauss->GaussPoint(ig); - gauss->SynchronizeGaussBase(gauss_tria); - - GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss); - tria->GetBSSA(&B[0][0], &xyz_list[0][0], gauss_tria); - tria->GetBprimeSSA(&Bprime[0][0], &xyz_list[0][0], gauss_tria); - - if(approximation==SSAHOApproximationEnum){ - this->StrainRateHO(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input); - this->StrainRateHO(&oldepsilon[0],&xyz_list[0][0],gauss,vxold_input,vyold_input); - material->GetViscosity3d(&viscosity, &epsilon[0]); - material->GetViscosity3d(&oldviscosity, &oldepsilon[0]); - - newviscosity=viscosity+viscosity_overshoot*(viscosity-oldviscosity); - } - else if (approximation==SSAFSApproximationEnum){ - this->StrainRateFS(&epsilons[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input); - material->GetViscosity3dFS(&newviscosity,&epsilons[0]); - } - else _error_("approximation " << approximation << " (" << EnumToStringx(approximation) << ") not supported yet"); - - D_scalar=2*newviscosity*gauss->weight*Jdet; - for (i=0;i<3;i++) D[i][i]=D_scalar; - - TripleMultiply( &B[0][0],3,numdof2d,1, - &D[0][0],3,3,0, - &Bprime[0][0],3,numdof2d,0, - &Ke_gg_gaussian[0][0],0); - - for(i=0;ivalues[i*numdof2d+j]+=Ke_gg_gaussian[i][j]; - } - - /*Transform Coordinate System*/ - ::TransformStiffnessMatrixCoord(Ke,tria->nodes,NUMVERTICES2D,XYEnum); - - /*Clean up and return*/ - delete tria->material; - delete tria; - delete gauss_tria; - delete gauss; - return Ke; -} -/*}}}*/ -/*FUNCTION Penta::CreateKMatrixStressbalanceSSA3dFriction{{{*/ -ElementMatrix* Penta::CreateKMatrixStressbalanceSSA3dFriction(void){ - - /*Initialize Element matrix and return if necessary*/ - if(IsFloating() || !IsOnBed()) return NULL; - - /*Build a tria element using the 3 nodes of the base of the penta. Then use - * the tria functionality to build a friction stiffness matrix on these 3 - * nodes: */ - Tria* tria=(Tria*)SpawnTria(0); //lower face is 0, upper face is 1. - ElementMatrix* Ke=tria->CreateKMatrixStressbalanceSSAFriction(); - delete tria->material; delete tria; - - /*clean-up and return*/ - return Ke; -} -/*}}}*/ -/*FUNCTION Penta::CreateKMatrixStressbalanceSSAHO{{{*/ -ElementMatrix* Penta::CreateKMatrixStressbalanceSSAHO(void){ - - /*compute all stiffness matrices for this element*/ - ElementMatrix* Ke1=CreateKMatrixStressbalanceSSA3d(); - ElementMatrix* Ke2=CreateKMatrixStressbalanceHO(); - ElementMatrix* Ke3=CreateKMatrixCouplingSSAHO(); - ElementMatrix* Ke =new ElementMatrix(Ke1,Ke2,Ke3); - - /*clean-up and return*/ - delete Ke1; - delete Ke2; - delete Ke3; - return Ke; -} -/*}}}*/ -/*FUNCTION Penta::CreateKMatrixStressbalanceSSAFS{{{*/ -ElementMatrix* Penta::CreateKMatrixStressbalanceSSAFS(void){ - - /*compute all stiffness matrices for this element*/ - ElementMatrix* Ke1=CreateKMatrixStressbalanceFS(); - int indices[3]={18,19,20}; - Ke1->StaticCondensation(3,&indices[0]); - int init = this->element_type; - this->element_type=P1Enum; //P1 needed for HO - ElementMatrix* Ke2=CreateKMatrixStressbalanceSSA3d(); - this->element_type=init; - ElementMatrix* Ke3=CreateKMatrixCouplingSSAFS(); - ElementMatrix* Ke =new ElementMatrix(Ke1,Ke2,Ke3); - - /*clean-up and return*/ - delete Ke1; - delete Ke2; - delete Ke3; - return Ke; -} -/*}}}*/ -/*FUNCTION Penta::CreateKMatrixStressbalanceL1L2{{{*/ -ElementMatrix* Penta::CreateKMatrixStressbalanceL1L2(void){ - - /*compute all stiffness matrices for this element*/ - ElementMatrix* Ke1=CreateKMatrixStressbalanceL1L2Viscous(); - ElementMatrix* Ke2=CreateKMatrixStressbalanceL1L2Friction(); - ElementMatrix* Ke =new ElementMatrix(Ke1,Ke2); - - /*clean-up and return*/ - delete Ke1; - delete Ke2; - return Ke; -} -/*}}}*/ -/*FUNCTION Penta::CreateKMatrixStressbalanceL1L2Viscous{{{*/ -ElementMatrix* Penta::CreateKMatrixStressbalanceL1L2Viscous(void){ - - /*Constants*/ - const int numdof2d=2*NUMVERTICES2D; - - /*Intermediaries */ - int i,j; - IssmDouble Jdet,viscosity; - IssmDouble xyz_list[NUMVERTICES][3]; - IssmDouble B[3][numdof2d]; - IssmDouble Bprime[3][numdof2d]; - IssmDouble Ke_gg_gaussian[numdof2d][numdof2d]; - IssmDouble D[3][3]= {0.0}; // material matrix, simple scalar matrix. - Tria *tria = NULL; - Penta *pentabase = NULL; - GaussPenta *gauss = NULL; - GaussTria *gauss_tria = NULL; - - /*Find penta on bed as this is a SSA elements: */ - pentabase=(Penta*)GetBasalElement(); - tria=pentabase->SpawnTria(0); //lower face is 0, upper face is 1. - - /*Initialize Element matrix*/ - ElementMatrix* Ke=new ElementMatrix(tria->nodes,NUMVERTICES2D,this->parameters,L1L2ApproximationEnum); - - /*Retrieve all inputs and parameters*/ - ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); - Input* vx_input=inputs->GetInput(VxEnum); _assert_(vx_input); - Input* vy_input=inputs->GetInput(VyEnum); _assert_(vy_input); - Input* surf_input=inputs->GetInput(SurfaceEnum); _assert_(surf_input); - - /* Start looping on the number of gaussian points: */ - gauss=new GaussPenta(5,5); - gauss_tria=new GaussTria(); - for(int ig=gauss->begin();igend();ig++){ - - gauss->GaussPoint(ig); - gauss->SynchronizeGaussBase(gauss_tria); - - GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss); - tria->GetBSSA(&B[0][0], &xyz_list[0][0], gauss_tria); - tria->GetBprimeSSA(&Bprime[0][0], &xyz_list[0][0], gauss_tria); - - /*Get viscosity for L1L2 model*/ - ViscosityL1L2(&viscosity,&xyz_list[0][0],gauss,vx_input,vy_input,surf_input); - - for(i=0;i<3;i++) D[i][i]=2*viscosity*gauss->weight*Jdet; - - TripleMultiply( &B[0][0],3,numdof2d,1, - &D[0][0],3,3,0, - &Bprime[0][0],3,numdof2d,0, - &Ke_gg_gaussian[0][0],0); - for(i=0;ivalues[i*numdof2d+j]+=Ke_gg_gaussian[i][j]; - } - - /*Transform Coordinate System*/ - ::TransformStiffnessMatrixCoord(Ke,tria->nodes,NUMVERTICES2D,XYEnum); - - /*Clean up and return*/ - delete tria->material; - delete tria; - delete gauss_tria; - delete gauss; - return Ke; -} -/*}}}*/ -/*FUNCTION Penta::CreateKMatrixStressbalanceL1L2Friction{{{*/ -ElementMatrix* Penta::CreateKMatrixStressbalanceL1L2Friction(void){ - - /*Initialize Element matrix and return if necessary*/ - if(IsFloating() || !IsOnBed()) return NULL; - - /*Build a tria element using the 3 nodes of the base of the penta. Then use - * the tria functionality to build a friction stiffness matrix on these 3 - * nodes: */ - Tria* tria=(Tria*)SpawnTria(0); //lower face is 0, upper face is 1. - ElementMatrix* Ke=tria->CreateKMatrixStressbalanceSSAFriction(); - delete tria->material; delete tria; - - /*clean-up and return*/ - return Ke; -} -/*}}}*/ -/*FUNCTION Penta::CreateKMatrixStressbalanceHO{{{*/ -ElementMatrix* Penta::CreateKMatrixStressbalanceHO(void){ - - /*compute all stiffness matrices for this element*/ - ElementMatrix* Ke1=CreateKMatrixStressbalanceHOViscous(); - ElementMatrix* Ke2=CreateKMatrixStressbalanceHOFriction(); - ElementMatrix* Ke =new ElementMatrix(Ke1,Ke2); - - /*clean-up and return*/ - delete Ke1; - delete Ke2; - return Ke; - -} -/*}}}*/ -/*FUNCTION Penta::CreateKMatrixStressbalanceHOViscous{{{*/ -ElementMatrix* Penta::CreateKMatrixStressbalanceHOViscous(void){ - - /*Intermediaries */ - int approximation; - IssmDouble xyz_list[NUMVERTICES][3]; - IssmDouble Jdet; - IssmDouble viscosity,oldviscosity,newviscosity,viscosity_overshoot; //viscosity - IssmDouble epsilon[5],oldepsilon[5]; /* epsilon=[exx,eyy,exy,exz,eyz];*/ - IssmDouble D_scalar; - GaussPenta *gauss=NULL; - - /*Fetch number of nodes and dof for this finite element*/ - int numnodes = this->NumberofNodes(); - int numdof = numnodes*NDOF2; - - /*Initialize Element matrix and vectors*/ - ElementMatrix* Ke = new ElementMatrix(nodes,numnodes,this->parameters,HOApproximationEnum); - IssmDouble* B = xNew(5*numdof); - IssmDouble* Bprime = xNew(5*numdof); - IssmDouble* D = xNewZeroInit(5*5); - - /*Retrieve all inputs and parameters*/ - inputs->GetInputValue(&approximation,ApproximationEnum); - ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); - this->parameters->FindParam(&viscosity_overshoot,StressbalanceViscosityOvershootEnum); - Input* vx_input=inputs->GetInput(VxEnum); _assert_(vx_input); - Input* vy_input=inputs->GetInput(VyEnum); _assert_(vy_input); - Input* vxold_input=inputs->GetInput(VxPicardEnum); _assert_(vxold_input); - Input* vyold_input=inputs->GetInput(VyPicardEnum); _assert_(vyold_input); - - /* Start looping on the number of gaussian points: */ - gauss=new GaussPenta(5,5); - for(int ig=gauss->begin();igend();ig++){ - - gauss->GaussPoint(ig); - - GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss); - GetBHO(&B[0], &xyz_list[0][0], gauss); - GetBprimeHO(&Bprime[0], &xyz_list[0][0], gauss); - - this->StrainRateHO(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input); - this->StrainRateHO(&oldepsilon[0],&xyz_list[0][0],gauss,vxold_input,vyold_input); - material->GetViscosity3d(&viscosity, &epsilon[0]); - material->GetViscosity3d(&oldviscosity, &oldepsilon[0]); - newviscosity=viscosity+viscosity_overshoot*(viscosity-oldviscosity); - - D_scalar=2*newviscosity*gauss->weight*Jdet; - for(int i=0;i<5;i++) D[i*5+i]=D_scalar; - - TripleMultiply(B,5,numdof,1, - D,5,5,0, - Bprime,5,numdof,0, - &Ke->values[0],1); - } - - /*Transform Coordinate System*/ - ::TransformStiffnessMatrixCoord(Ke,nodes,numnodes,XYEnum); - - /*Clean up and return*/ - delete gauss; - xDelete(D); - xDelete(Bprime); - xDelete(B); - return Ke; -} -/*}}}*/ -/*FUNCTION Penta::CreateKMatrixStressbalanceHOFriction{{{*/ -ElementMatrix* Penta::CreateKMatrixStressbalanceHOFriction(void){ - - /*Intermediaries */ - bool mainlyfloating; - int i,j; - int analysis_type,migration_style; - int point1; - IssmDouble xyz_list[NUMVERTICES][3]; - IssmDouble xyz_list_tria[NUMVERTICES2D][3]; - IssmDouble alpha2,Jdet; - IssmDouble fraction1,fraction2; - IssmDouble gllevelset; - IssmDouble phi=1.0; - IssmDouble DL_scalar; - Friction *friction = NULL; - GaussPenta *gauss = NULL; - - /*Initialize Element matrix and return if necessary*/ - if(IsFloating() || !IsOnBed()) return NULL; - - /*Fetch number of nodes and dof for this finite element*/ - int numnodes = this->NumberofNodes(); - int numdof = numnodes*NDOF2; - - /*Initialize Element matrix and vectors*/ - ElementMatrix* Ke = new ElementMatrix(nodes,numnodes,this->parameters,HOApproximationEnum); - IssmDouble* B = xNew(2*numdof); - IssmDouble* D = xNewZeroInit(2*2); - - /*Retrieve all inputs and parameters*/ - ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); - for(i=0;iFindParam(&migration_style,GroundinglineMigrationEnum); - parameters->FindParam(&analysis_type,AnalysisTypeEnum); - Input* surface_input=inputs->GetInput(SurfaceEnum); _assert_(surface_input); - Input* vx_input=inputs->GetInput(VxEnum); _assert_(vx_input); - Input* vy_input=inputs->GetInput(VyEnum); _assert_(vy_input); - Input* vz_input=inputs->GetInput(VzEnum); _assert_(vz_input); - Input* gllevelset_input=NULL; - - /*build friction object, used later on: */ - friction=new Friction(this,2); - - /*Recover portion of element that is grounded*/ - if(migration_style==SubelementMigrationEnum) phi=this->GetGroundedPortion(&xyz_list_tria[0][0]); - if(migration_style==SubelementMigration2Enum){ - gllevelset_input=inputs->GetInput(MaskGroundediceLevelsetEnum); _assert_(gllevelset_input); - this->GetGroundedPart(&point1,&fraction1,&fraction2,&mainlyfloating); - //gauss=new GaussPenta(point1,fraction1,fraction2,mainlyfloating,2); - gauss=new GaussPenta(0,1,2,2); - } - else{ - gauss=new GaussPenta(0,1,2,2); - } - - /* Start looping on the number of gaussian points: */ - for(int ig=gauss->begin();igend();ig++){ - - gauss->GaussPoint(ig); - - GetTriaJacobianDeterminant(&Jdet,&xyz_list_tria[0][0],gauss); - GetBHOFriction(&B[0],gauss); - - friction->GetAlpha2(&alpha2,gauss,vx_input,vy_input,vz_input); - if(migration_style==SubelementMigrationEnum) alpha2=phi*alpha2; - if(migration_style==SubelementMigration2Enum){ - gllevelset_input->GetInputValue(&gllevelset, gauss); - if(gllevelset<0) alpha2=0; - } - - DL_scalar=alpha2*gauss->weight*Jdet; - for (i=0;i<2;i++) D[i*2+i]=DL_scalar; - - TripleMultiply(B,2,numdof,1, - D,2,2,0, - B,2,numdof,0, - &Ke->values[0],1); - } - - /*Transform Coordinate System*/ - ::TransformStiffnessMatrixCoord(Ke,nodes,numnodes,XYEnum); - - /*Clean up and return*/ - delete gauss; - xDelete(D); - xDelete(B); - delete friction; - return Ke; -} -/*}}}*/ -/*FUNCTION Penta::CreateKMatrixStressbalanceHOFS{{{*/ -ElementMatrix* Penta::CreateKMatrixStressbalanceHOFS(void){ - - /*compute all stiffness matrices for this element*/ - ElementMatrix* Ke1=CreateKMatrixStressbalanceFS(); - int indices[3]={18,19,20}; - Ke1->StaticCondensation(3,&indices[0]); - int init = this->element_type; - this->element_type=P1Enum; //P1 needed for HO - ElementMatrix* Ke2=CreateKMatrixStressbalanceHO(); - this->element_type=init; - ElementMatrix* Ke3=CreateKMatrixCouplingHOFS(); - ElementMatrix* Ke =new ElementMatrix(Ke1,Ke2,Ke3); - - /*clean-up and return*/ - delete Ke1; - delete Ke2; - delete Ke3; - return Ke; -} -/*}}}*/ -/*FUNCTION Penta::CreateKMatrixStressbalanceFS{{{*/ -ElementMatrix* Penta::CreateKMatrixStressbalanceFS(void){ - - ElementMatrix* Ke1 = NULL; - ElementMatrix* Ke2 = NULL; - ElementMatrix* Ke = NULL; - - /*compute all stiffness matrices for this element*/ - Ke1=CreateKMatrixStressbalanceFSViscous(); - Ke2=CreateKMatrixStressbalanceFSFriction(); - Ke =new ElementMatrix(Ke1,Ke2); - - /*clean-up and return*/ - delete Ke1; - delete Ke2; - return Ke; - -} -/*}}}*/ /*FUNCTION Penta::KMatrixGLSstabilization{{{*/ void Penta::KMatrixGLSstabilization(ElementMatrix* Ke){ @@ -6556,511 +4625,7 @@ } /*}}}*/ -/*FUNCTION Penta::CreateKMatrixStressbalanceFSViscous {{{*/ -ElementMatrix* Penta::CreateKMatrixStressbalanceFSViscous(void){ - - /*Intermediaries */ - int i,approximation; - IssmDouble Jdet,viscosity,FSreconditioning,D_scalar; - IssmDouble xyz_list[NUMVERTICES][3]; - IssmDouble epsilon[6]; /* epsilon=[exx,eyy,ezz,exy,exz,eyz];*/ - GaussPenta *gauss=NULL; - - /*If on water or not FS, skip stiffness: */ - inputs->GetInputValue(&approximation,ApproximationEnum); - if(approximation!=FSApproximationEnum && approximation!=SSAFSApproximationEnum && approximation!=HOFSApproximationEnum) return NULL; - - /*Fetch number of nodes and dof for this finite element*/ - int vnumnodes = this->NumberofNodesVelocity(); - int pnumnodes = this->NumberofNodesPressure(); - int numdof = vnumnodes*NDOF3 + pnumnodes*NDOF1; - - /*Prepare coordinate system list*/ - int* cs_list = xNew(vnumnodes+pnumnodes); - for(i=0;iparameters,FSvelocityEnum); - IssmDouble* B = xNew(8*numdof); - IssmDouble* Bprime = xNew(8*numdof); - IssmDouble* D = xNewZeroInit(8*8); - - /*Retrieve all inputs and parameters*/ - ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); - parameters->FindParam(&FSreconditioning,StressbalanceFSreconditioningEnum); - Input* vx_input=inputs->GetInput(VxEnum); _assert_(vx_input); - Input* vy_input=inputs->GetInput(VyEnum); _assert_(vy_input); - Input* vz_input=inputs->GetInput(VzEnum); _assert_(vz_input); - - /* Start looping on the number of gaussian points: */ - gauss=new GaussPenta(5,5); - for(int ig=gauss->begin();igend();ig++){ - - gauss->GaussPoint(ig); - - GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss); - GetBFS(B,&xyz_list[0][0],gauss); - GetBprimeFS(Bprime,&xyz_list[0][0],gauss); - - this->StrainRateFS(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input); - material->GetViscosity3dFS(&viscosity,&epsilon[0]); - - D_scalar=gauss->weight*Jdet; - for(i=0;i<6;i++) D[i*8+i] = +D_scalar*2.*viscosity; - for(i=6;i<8;i++) D[i*8+i] = -D_scalar*FSreconditioning; - - TripleMultiply(B,8,numdof,1, - D,8,8,0, - Bprime,8,numdof,0, - Ke->values,1); - } - - /*Transform Coordinate System*/ - ::TransformStiffnessMatrixCoord(Ke,nodes,(vnumnodes+pnumnodes),cs_list); - - /*Clean up and return*/ - delete gauss; - xDelete(cs_list); - xDelete(B); - xDelete(Bprime); - xDelete(D); - return Ke; -} -/*}}}*/ -/*FUNCTION Penta::CreateKMatrixStressbalanceFSFriction{{{*/ -ElementMatrix* Penta::CreateKMatrixStressbalanceFSFriction(void){ - - /*Intermediaries */ - int i,j; - int analysis_type,approximation; - IssmDouble alpha2,Jdet2d; - IssmDouble xyz_list[NUMVERTICES][3]; - IssmDouble xyz_list_tria[NUMVERTICES2D][3]; - Friction *friction = NULL; - GaussPenta *gauss = NULL; - - /*Initialize Element matrix and return if necessary*/ - if(IsFloating() || !IsOnBed()) return NULL; - - /*If on water or not FS, skip stiffness: */ - inputs->GetInputValue(&approximation,ApproximationEnum); - if(approximation!=FSApproximationEnum && approximation!=SSAFSApproximationEnum && approximation!=HOFSApproximationEnum) return NULL; - - /*Fetch number of nodes and dof for this finite element*/ - int vnumnodes = this->NumberofNodesVelocity(); - int pnumnodes = this->NumberofNodesPressure(); - int numdof = vnumnodes*NDOF3 + pnumnodes*NDOF1; - - /*Initialize Element matrix and vectors*/ - ElementMatrix* Ke = new ElementMatrix(nodes,vnumnodes+pnumnodes,this->parameters,FSvelocityEnum); - IssmDouble* BFriction = xNew(2*numdof); - IssmDouble* D = xNewZeroInit(2*2); - - /*Retrieve all inputs and parameters*/ - ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); - parameters->FindParam(&analysis_type,AnalysisTypeEnum); - Input* vx_input=inputs->GetInput(VxEnum); _assert_(vx_input); - Input* vy_input=inputs->GetInput(VyEnum); _assert_(vy_input); - Input* vz_input=inputs->GetInput(VzEnum); _assert_(vz_input); - for(i=0;ibegin();igend();ig++){ - - gauss->GaussPoint(ig); - - GetTriaJacobianDeterminant(&Jdet2d, &xyz_list_tria[0][0],gauss); - GetLFS(BFriction,gauss); - - friction->GetAlpha2(&alpha2,gauss,vx_input,vy_input,vz_input); - - D[0*2+0] = +alpha2*gauss->weight*Jdet2d; //taub_x = -alpha2 vx - D[1*2+1] = +alpha2*gauss->weight*Jdet2d; //taub_y = -alpha2 vy - - TripleMultiply(BFriction,2,numdof,1, - D,2,2,0, - BFriction,2,numdof,0, - Ke->values,1); - } - - /*DO NOT Transform Coordinate System: this stiffness matrix is already expressed in tangential coordinates*/ - //::TransformStiffnessMatrixCoord(Ke,nodes,NUMVERTICES,XYZEnum); - - /*Clean up and return*/ - delete gauss; - delete friction; - xDelete(BFriction); - xDelete(D); - return Ke; -} -/*}}}*/ -/*FUNCTION Penta::CreateKMatrixStressbalanceVert {{{*/ -ElementMatrix* Penta::CreateKMatrixStressbalanceVert(void){ - - /*compute all stiffness matrices for this element*/ - ElementMatrix* Ke1=CreateKMatrixStressbalanceVertVolume(); - ElementMatrix* Ke2=CreateKMatrixStressbalanceVertSurface(); - ElementMatrix* Ke =new ElementMatrix(Ke1,Ke2); - - /*clean-up and return*/ - delete Ke1; - delete Ke2; - return Ke; - -} -/*}}}*/ -/*FUNCTION Penta::CreateKMatrixStressbalanceVertVolume {{{*/ -ElementMatrix* Penta::CreateKMatrixStressbalanceVertVolume(void){ - - /*Intermediaries */ - IssmDouble Jdet; - IssmDouble xyz_list[NUMVERTICES][3]; - IssmDouble B[NDOF1][NUMVERTICES]; - IssmDouble Bprime[NDOF1][NUMVERTICES]; - IssmDouble DL_scalar; - GaussPenta *gauss=NULL; - - /*Initialize Element matrix*/ - ElementMatrix* Ke=new ElementMatrix(nodes,NUMVERTICES,this->parameters,NoneApproximationEnum); - - /*Retrieve all inputs and parameters*/ - ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); - - /* Start looping on the number of gaussian points: */ - gauss=new GaussPenta(2,2); - for(int ig=gauss->begin();igend();ig++){ - - gauss->GaussPoint(ig); - - GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss); - GetBVert(&B[0][0], &xyz_list[0][0], gauss); - GetBprimeVert(&Bprime[0][0], &xyz_list[0][0], gauss); - - DL_scalar=gauss->weight*Jdet; - - TripleMultiply( &B[0][0],1,NUMVERTICES,1, - &DL_scalar,1,1,0, - &Bprime[0][0],1,NUMVERTICES,0, - &Ke->values[0],1); - } - - /*Clean up and return*/ - delete gauss; - return Ke; -} -/*}}}*/ -/*FUNCTION Penta::CreateKMatrixStressbalanceVertSurface {{{*/ -ElementMatrix* Penta::CreateKMatrixStressbalanceVertSurface(void){ - - if (!IsOnSurface()) return NULL; - - /*Constants*/ - const int numdof=NDOF1*NUMVERTICES; - - /*Intermediaries */ - int i,j; - IssmDouble xyz_list[NUMVERTICES][3]; - IssmDouble xyz_list_tria[NUMVERTICES2D][3]; - IssmDouble surface_normal[3]; - IssmDouble Jdet2d,DL_scalar; - IssmDouble basis[NUMVERTICES]; - GaussPenta *gauss=NULL; - - /*Initialize Element matrix*/ - ElementMatrix* Ke=new ElementMatrix(nodes,NUMVERTICES,this->parameters,NoneApproximationEnum); - - /*Retrieve all inputs and parameters*/ - ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); - for(i=0;ibegin();igend();ig++){ - - gauss->GaussPoint(ig); - - GetTriaJacobianDeterminant(&Jdet2d, &xyz_list_tria[0][0],gauss); - GetNodalFunctionsP1(&basis[0], gauss); - - DL_scalar= - gauss->weight*Jdet2d*surface_normal[2]; - - TripleMultiply( basis,1,numdof,1, - &DL_scalar,1,1,0, - basis,1,numdof,0, - &Ke->values[0],1); - } - - /*Clean up and return*/ - delete gauss; - return Ke; -} -/*}}}*/ -/*FUNCTION Penta::CreateJacobianStressbalanceHoriz{{{*/ -ElementMatrix* Penta::CreateJacobianStressbalanceHoriz(void){ - - int approximation; - inputs->GetInputValue(&approximation,ApproximationEnum); - - switch(approximation){ - case SSAApproximationEnum: - return CreateJacobianStressbalanceSSA2d(); - case HOApproximationEnum: - return CreateJacobianStressbalanceHO(); - case FSApproximationEnum: - return CreateJacobianStressbalanceFS(); - case NoneApproximationEnum: - return NULL; - default: - _error_("Approximation " << EnumToStringx(approximation) << " not supported yet"); - } -} -/*}}}*/ -/*FUNCTION Penta::CreateJacobianStressbalanceSSA2d{{{*/ -ElementMatrix* Penta::CreateJacobianStressbalanceSSA2d(void){ - - /*Figure out if this penta is collapsed. If so, then bailout, except if it is at the - bedrock, in which case we spawn a tria element using the 3 first nodes, and use it to build - the stiffness matrix. */ - if (!IsOnBed()) return NULL; - - /*Depth average some fields*/ - switch(this->material->ObjectEnum()){ - case MaticeEnum: - this->InputDepthAverageAtBase(MaterialsRheologyBEnum,MaterialsRheologyBbarEnum,MaterialsEnum); - this->InputDepthAverageAtBase(DamageDEnum,DamageDbarEnum,MaterialsEnum); - break; - default: - _error_("material "<material->ObjectEnum())<<" not supported"); - } - - /*Call Tria function*/ - Tria* tria=(Tria*)SpawnTria(0); //lower face is 0, upper face is 1. - ElementMatrix* Ke=tria->CreateJacobianStressbalanceSSA(); - delete tria->material; delete tria; - - /*Delete averaged inputs*/ - this->material->inputs->DeleteInput(MaterialsRheologyBbarEnum); - this->material->inputs->DeleteInput(DamageDbarEnum); - - /*clean up and return*/ - return Ke; -} -/*}}}*/ -/*FUNCTION Penta::CreateJacobianStressbalanceHO{{{*/ -ElementMatrix* Penta::CreateJacobianStressbalanceHO(void){ - - /*Intermediaries */ - int i,j; - IssmDouble xyz_list[NUMVERTICES][3]; - IssmDouble Jdet; - IssmDouble eps1dotdphii,eps1dotdphij; - IssmDouble eps2dotdphii,eps2dotdphij; - IssmDouble mu_prime; - IssmDouble epsilon[5]; /* epsilon=[exx,eyy,exy,exz,eyz];*/ - IssmDouble eps1[3],eps2[3]; - IssmDouble dphi[3][NUMVERTICES]; - GaussPenta *gauss=NULL; - - /*Initialize Jacobian with regular HO (first part of the Gateau derivative)*/ - ElementMatrix* Ke=CreateKMatrixStressbalanceHO(); - - /*Retrieve all inputs and parameters*/ - ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); - Input* vx_input=inputs->GetInput(VxEnum); _assert_(vx_input); - Input* vy_input=inputs->GetInput(VyEnum); _assert_(vy_input); - - /* Start looping on the number of gaussian points: */ - gauss=new GaussPenta(5,5); - for(int ig=gauss->begin();igend();ig++){ - - gauss->GaussPoint(ig); - - GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss); - GetNodalFunctionsP1Derivatives(&dphi[0][0],&xyz_list[0][0],gauss); - - this->StrainRateHO(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input); - material->GetViscosityDerivativeEpsSquare(&mu_prime,&epsilon[0]); - eps1[0]=2*epsilon[0]+epsilon[1]; eps2[0]=epsilon[2]; - eps1[1]=epsilon[2]; eps2[1]=epsilon[0]+2*epsilon[1]; - eps1[2]=epsilon[3]; eps2[2]=epsilon[4]; - - for(i=0;i<6;i++){ - for(j=0;j<6;j++){ - eps1dotdphii=eps1[0]*dphi[0][i]+eps1[1]*dphi[1][i]+eps1[2]*dphi[2][i]; - eps1dotdphij=eps1[0]*dphi[0][j]+eps1[1]*dphi[1][j]+eps1[2]*dphi[2][j]; - eps2dotdphii=eps2[0]*dphi[0][i]+eps2[1]*dphi[1][i]+eps2[2]*dphi[2][i]; - eps2dotdphij=eps2[0]*dphi[0][j]+eps2[1]*dphi[1][j]+eps2[2]*dphi[2][j]; - - Ke->values[12*(2*i+0)+2*j+0]+=gauss->weight*Jdet*2*mu_prime*eps1dotdphij*eps1dotdphii; - Ke->values[12*(2*i+0)+2*j+1]+=gauss->weight*Jdet*2*mu_prime*eps2dotdphij*eps1dotdphii; - Ke->values[12*(2*i+1)+2*j+0]+=gauss->weight*Jdet*2*mu_prime*eps1dotdphij*eps2dotdphii; - Ke->values[12*(2*i+1)+2*j+1]+=gauss->weight*Jdet*2*mu_prime*eps2dotdphij*eps2dotdphii; - } - } - } - - /*Transform Coordinate System*/ - ::TransformStiffnessMatrixCoord(Ke,nodes,NUMVERTICES,XYEnum); - - /*Clean up and return*/ - delete gauss; - return Ke; -} -/*}}}*/ -/*FUNCTION Penta::CreateJacobianStressbalanceFS{{{*/ -ElementMatrix* Penta::CreateJacobianStressbalanceFS(void){ - - /*Intermediaries */ - int i,j; - IssmDouble xyz_list[NUMVERTICES][3]; - IssmDouble Jdet; - IssmDouble eps1dotdphii,eps1dotdphij; - IssmDouble eps2dotdphii,eps2dotdphij; - IssmDouble eps3dotdphii,eps3dotdphij; - IssmDouble mu_prime; - IssmDouble epsilon[6]; /* epsilon=[exx,eyy,ezz,exy,exz,eyz];*/ - IssmDouble eps1[3],eps2[3],eps3[3]; - GaussPenta *gauss=NULL; - - /*If on water or not FS, skip stiffness: */ - //inputs->GetInputValue(&approximation,ApproximationEnum); - //if(approximation!=FSApproximationEnum && approximation!=SSAFSApproximationEnum && approximation!=HOFSApproximationEnum) return NULL; - - /*Fetch number of nodes and dof for this finite element*/ - int vnumnodes = this->NumberofNodesVelocity(); - int pnumnodes = this->NumberofNodesPressure(); - int numdof = vnumnodes*NDOF3 + pnumnodes*NDOF1; - - /*Prepare coordinate system list*/ - int* cs_list = xNew(vnumnodes+pnumnodes); - for(i=0;i(3*vnumnodes); - - /*Retrieve all inputs and parameters*/ - ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); - Input* vx_input=inputs->GetInput(VxEnum); _assert_(vx_input); - Input* vy_input=inputs->GetInput(VyEnum); _assert_(vy_input); - Input* vz_input=inputs->GetInput(VzEnum); _assert_(vz_input); - - /* Start looping on the number of gaussian points: */ - gauss=new GaussPenta(5,5); - for(int ig=gauss->begin();igend();ig++){ - - gauss->GaussPoint(ig); - - GetJacobianDeterminant(&Jdet,&xyz_list[0][0],gauss); - GetNodalFunctionsDerivativesVelocity(dbasis,&xyz_list[0][0],gauss); - - //this->StrainRateFS(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input); - //material->GetViscosityDerivativeEpsSquareFS(&mu_prime,&epsilon[0]); - //eps1[0]=epsilon[0]; eps2[0]=epsilon[3]; eps3[0]=epsilon[4]; - //eps1[1]=epsilon[3]; eps2[1]=epsilon[1]; eps3[1]=epsilon[5]; - //eps1[2]=epsilon[4]; eps2[2]=epsilon[5]; eps3[2]=epsilon[2]; - this->StrainRateHO(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input); - material->GetViscosityDerivativeEpsSquare(&mu_prime,&epsilon[0]); - eps1[0]=epsilon[0]; eps2[0]=epsilon[2]; eps3[0]=epsilon[3]; - eps1[1]=epsilon[2]; eps2[1]=epsilon[1]; eps3[1]=epsilon[4]; - eps1[2]=epsilon[3]; eps2[2]=epsilon[4]; eps3[2]= -epsilon[0] -epsilon[1]; - - for(i=0;ivalues[numdof*(3*i+0)+3*j+0]+=gauss->weight*Jdet*2*mu_prime*eps1dotdphij*eps1dotdphii; - Ke->values[numdof*(3*i+0)+3*j+1]+=gauss->weight*Jdet*2*mu_prime*eps2dotdphij*eps1dotdphii; - Ke->values[numdof*(3*i+0)+3*j+2]+=gauss->weight*Jdet*2*mu_prime*eps3dotdphij*eps1dotdphii; - - Ke->values[numdof*(3*i+1)+3*j+0]+=gauss->weight*Jdet*2*mu_prime*eps1dotdphij*eps2dotdphii; - Ke->values[numdof*(3*i+1)+3*j+1]+=gauss->weight*Jdet*2*mu_prime*eps2dotdphij*eps2dotdphii; - Ke->values[numdof*(3*i+1)+3*j+2]+=gauss->weight*Jdet*2*mu_prime*eps3dotdphij*eps2dotdphii; - - Ke->values[numdof*(3*i+2)+3*j+0]+=gauss->weight*Jdet*2*mu_prime*eps1dotdphij*eps3dotdphii; - Ke->values[numdof*(3*i+2)+3*j+1]+=gauss->weight*Jdet*2*mu_prime*eps2dotdphij*eps3dotdphii; - Ke->values[numdof*(3*i+2)+3*j+2]+=gauss->weight*Jdet*2*mu_prime*eps3dotdphij*eps3dotdphii; - } - } - } - - /*Transform Coordinate System*/ - ::TransformStiffnessMatrixCoord(Ke,nodes,(vnumnodes+pnumnodes),cs_list); - - /*Clean up and return*/ - xDelete(cs_list); - xDelete(dbasis); - delete gauss; - return Ke; -} -/*}}}*/ #endif -#ifdef _HAVE_BALANCED_ -/*FUNCTION Penta::CreateKMatrixBalancethickness {{{*/ -ElementMatrix* Penta::CreateKMatrixBalancethickness(void){ - - /*Figure out if this penta is collapsed. If so, then bailout, except if it is at the - bedrock, in which case we spawn a tria element using the 3 first nodes, and use it to build - the stiffness matrix. */ - if (!IsOnBed()) return NULL; - - /*Depth Averaging Vx and Vy*/ - this->InputDepthAverageAtBase(VxEnum,VxAverageEnum); - this->InputDepthAverageAtBase(VyEnum,VyAverageEnum); - - /*Spawn Tria element from the base of the Penta: */ - Tria* tria=(Tria*)SpawnTria(0); //lower face is 0, upper face is 1. - ElementMatrix* Ke=tria->CreateKMatrixBalancethickness(); - delete tria->material; delete tria; - - /*Delete Vx and Vy averaged*/ - this->inputs->DeleteInput(VxAverageEnum); - this->inputs->DeleteInput(VyAverageEnum); - - /*clean up and return*/ - return Ke; -} -/*}}}*/ -#endif - #ifdef _HAVE_HYDROLOGY_ -/*FUNCTION Penta::CreateKMatrixHydrologyDCInefficient {{{*/ -ElementMatrix* Penta::CreateKMatrixHydrologyDCInefficient(void){ - - if (!IsOnBed()) return NULL; - - Tria* tria=(Tria*)SpawnTria(0); //lower face is 0, upper face is 1. - ElementMatrix* Ke=tria->CreateKMatrixHydrologyDCInefficient(); - delete tria->material; delete tria; - - /*clean up and return*/ - return Ke; -} -/*}}}*/ -/*FUNCTION Penta::CreateKMatrixHydrologyDCEfficient {{{*/ -ElementMatrix* Penta::CreateKMatrixHydrologyDCEfficient(void){ - - if (!IsOnBed()) return NULL; - - Tria* tria=(Tria*)SpawnTria(0); //lower face is 0, upper face is 1. - ElementMatrix* Ke=tria->CreateKMatrixHydrologyDCEfficient(); - delete tria->material; delete tria; - - /*clean up and return*/ - return Ke; -} -/*}}}*/ /*FUNCTION Penta::CreateEPLDomainMassMatrix {{{*/ ElementMatrix* Penta::CreateEPLDomainMassMatrix(void){ Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Penta.h (revision 16992) +++ /issm/trunk-jpl/src/c/classes/Elements/Penta.h (revision 16993) @@ -72,5 +72,4 @@ void SetwiseNodeConnectivity(int* d_nz,int* o_nz,Node* node,bool* flags,int* flagsindices,int set1_enum,int set2_enum); void CreateDVector(Vector* df); - void CreateJacobianMatrix(Matrix* Jff); void Delta18oParameterization(void); void EnthalpyToThermal(IssmDouble* ptemperature,IssmDouble* pwaterfraction,IssmDouble enthalpy,IssmDouble pressure); @@ -195,8 +194,4 @@ void NormalTop(IssmDouble* bed_normal, IssmDouble* xyz_list); ElementMatrix* CreateBasalMassMatrix(void); - ElementMatrix* CreateKMatrix(void); - ElementMatrix* CreateKMatrixMasstransport(void); - ElementMatrix* CreateKMatrixFreeSurfaceTop(void); - ElementMatrix* CreateKMatrixFreeSurfaceBase(void); IssmDouble EnthalpyDiffusionParameter(IssmDouble enthalpy,IssmDouble pressure); IssmDouble EnthalpyDiffusionParameterVolume(int numvertices,IssmDouble* enthalpy,IssmDouble* pressure); @@ -260,51 +255,10 @@ #ifdef _HAVE_STRESSBALANCE_ - ElementMatrix* CreateKMatrixCouplingSSAHO(void); - ElementMatrix* CreateKMatrixCouplingSSAHOViscous(void); - ElementMatrix* CreateKMatrixCouplingSSAHOFriction(void); - ElementMatrix* CreateKMatrixCouplingSSAFS(void); - ElementMatrix* CreateKMatrixCouplingSSAFSViscous(void); - ElementMatrix* CreateKMatrixCouplingSSAFSFriction(void); - ElementMatrix* CreateKMatrixCouplingHOFS(void); - ElementMatrix* CreateKMatrixStressbalanceHoriz(void); - ElementMatrix* CreateKMatrixAdjointHoriz(void); ElementVector* CreateDVectorStressbalanceHoriz(void); ElementVector* CreateDVectorStressbalanceFS(void); - ElementMatrix* CreateKMatrixStressbalanceSIA(void); - ElementMatrix* CreateKMatrixStressbalanceSSA2d(void); - ElementMatrix* CreateKMatrixStressbalanceSSA3d(void); - ElementMatrix* CreateKMatrixStressbalanceSSA3dViscous(void); - ElementMatrix* CreateKMatrixStressbalanceSSA3dFriction(void); - ElementMatrix* CreateKMatrixStressbalanceSSAHO(void); - ElementMatrix* CreateKMatrixStressbalanceSSAFS(void); - ElementMatrix* CreateKMatrixStressbalanceL1L2(void); - ElementMatrix* CreateKMatrixStressbalanceL1L2Viscous(void); - ElementMatrix* CreateKMatrixStressbalanceL1L2Friction(void); - ElementMatrix* CreateKMatrixStressbalanceHO(void); - ElementMatrix* CreateKMatrixStressbalanceHOViscous(void); - ElementMatrix* CreateKMatrixStressbalanceHOFriction(void); - ElementMatrix* CreateKMatrixStressbalanceHOFS(void); - ElementMatrix* CreateKMatrixStressbalanceFS(void); - ElementMatrix* CreateKMatrixStressbalanceFSViscous(void); void KMatrixGLSstabilization(ElementMatrix* Ke); - ElementMatrix* CreateKMatrixStressbalanceFSFriction(void); - ElementMatrix* CreateKMatrixStressbalanceVert(void); - ElementMatrix* CreateKMatrixStressbalanceVertVolume(void); - ElementMatrix* CreateKMatrixStressbalanceVertSurface(void); - ElementMatrix* CreateJacobianStressbalanceHoriz(void); - ElementMatrix* CreateJacobianStressbalanceSSA2d(void); - ElementMatrix* CreateJacobianStressbalanceHO(void); - ElementMatrix* CreateJacobianStressbalanceFS(void); - #endif - - #ifdef _HAVE_CONTROL_ - ElementMatrix* CreateKMatrixAdjointSSA2d(void); - ElementMatrix* CreateKMatrixAdjointHO(void); - ElementMatrix* CreateKMatrixAdjointFS(void); #endif #ifdef _HAVE_HYDROLOGY_ - ElementMatrix* CreateKMatrixHydrologyDCInefficient(void); - ElementMatrix* CreateKMatrixHydrologyDCEfficient(void); ElementMatrix* CreateEPLDomainMassMatrix(void); void GetHydrologyDCInefficientHmax(IssmDouble* ph_max, Node* innode); @@ -319,19 +273,8 @@ void UpdateConstraintsExtrudeFromTop(void){_error_("not implemented yet");}; #ifdef _HAVE_THERMAL_ - ElementMatrix* CreateKMatrixEnthalpy(void); - ElementMatrix* CreateKMatrixEnthalpyVolume(void); - ElementMatrix* CreateKMatrixEnthalpyShelf(void); - ElementMatrix* CreateKMatrixThermal(void); - ElementMatrix* CreateKMatrixMelting(void); - ElementMatrix* CreateKMatrixThermalVolume(void); - ElementMatrix* CreateKMatrixThermalShelf(void); void UpdateBasalConstraintsEnthalpy(void); void ComputeBasalMeltingrate(void); void DrainWaterfraction(IssmDouble* drainrate_element); #endif - - #ifdef _HAVE_BALANCED_ - ElementMatrix* CreateKMatrixBalancethickness(void); - #endif /*}}}*/ }; Index: /issm/trunk-jpl/src/c/classes/Elements/Seg.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Seg.cpp (revision 16992) +++ /issm/trunk-jpl/src/c/classes/Elements/Seg.cpp (revision 16993) @@ -138,199 +138,4 @@ } /*}}}*/ -/*FUNCTION Seg::CreateKMatrixFreeSurfaceTop {{{*/ -ElementMatrix* Seg::CreateKMatrixFreeSurfaceTop(void){ - - /*Intermediaries */ - int stabilization; - IssmDouble Jdet,D_scalar,dt,h; - IssmDouble vx,vel; - IssmDouble xyz_list[NUMVERTICES][3]; - - /*Fetch number of nodes for this finite element*/ - int numnodes = this->NumberofNodes(); - - /*Initialize Element matrix and vectors*/ - ElementMatrix* Ke = new ElementMatrix(nodes,numnodes,this->parameters,NoneApproximationEnum); - IssmDouble* basis = xNew(numnodes); - IssmDouble* B = xNew(1*numnodes); - IssmDouble* Bprime = xNew(1*numnodes); - - /*Retrieve all inputs and parameters*/ - ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); - this->parameters->FindParam(&dt,TimesteppingTimeStepEnum); - this->parameters->FindParam(&stabilization,MasstransportStabilizationEnum); - Input* vx_input=inputs->GetInput(VxEnum); _assert_(vx_input); - h=this->GetSize(); - - /* Start looping on the number of gaussian points: */ - GaussSeg *gauss=new GaussSeg(2); - for(int ig=gauss->begin();igend();ig++){ - - gauss->GaussPoint(ig); - - GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss); - GetNodalFunctions(basis,gauss); - - vx_input->GetInputValue(&vx,gauss); - - D_scalar=gauss->weight*Jdet; - - TripleMultiply(basis,1,numnodes,1, - &D_scalar,1,1,0, - basis,1,numnodes,0, - &Ke->values[0],1); - - GetNodalFunctions(B,gauss); - GetBprimeMasstransport(Bprime,&xyz_list[0][0],gauss); - - D_scalar=dt*gauss->weight*Jdet*vx; - TripleMultiply(B,1,numnodes,1, - &D_scalar,1,1,0, - Bprime,1,numnodes,0, - &Ke->values[0],1); - - if(stabilization==2){ - /*Streamline upwinding*/ - vel=fabs(vx)+1.e-8; - D_scalar=dt*gauss->weight*Jdet*h/(2.*vel)*vx; - } - else if(stabilization==1){ - /*SSA*/ - vx_input->GetInputAverage(&vx); - D_scalar=dt*gauss->weight*Jdet*h/2.*fabs(vx); - } - if(stabilization==1 || stabilization==2){ - TripleMultiply(Bprime,1,numnodes,1, - &D_scalar,1,1,0, - Bprime,1,numnodes,0, - &Ke->values[0],1); - } - } - - /*Clean up and return*/ - xDelete(basis); - xDelete(B); - xDelete(Bprime); - delete gauss; - return Ke; -} -/*}}}*/ -/*FUNCTION Seg::CreateKMatrixFreeSurfaceBase {{{*/ -ElementMatrix* Seg::CreateKMatrixFreeSurfaceBase(void){ - - /*Intermediaries */ - int stabilization; - IssmDouble Jdet,D_scalar,dt,h; - IssmDouble vx,vel; - IssmDouble xyz_list[NUMVERTICES][3]; - - /*Fetch number of nodes for this finite element*/ - int numnodes = this->NumberofNodes(); - - /*Initialize Element matrix and vectors*/ - ElementMatrix* Ke = new ElementMatrix(nodes,numnodes,this->parameters,NoneApproximationEnum); - IssmDouble* basis = xNew(numnodes); - IssmDouble* B = xNew(1*numnodes); - IssmDouble* Bprime = xNew(1*numnodes); - - /*Retrieve all inputs and parameters*/ - ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); - this->parameters->FindParam(&dt,TimesteppingTimeStepEnum); - this->parameters->FindParam(&stabilization,MasstransportStabilizationEnum); - Input* vx_input=inputs->GetInput(VxEnum); _assert_(vx_input); - h=this->GetSize(); - - /* Start looping on the number of gaussian points: */ - GaussSeg *gauss=new GaussSeg(2); - for(int ig=gauss->begin();igend();ig++){ - - gauss->GaussPoint(ig); - - GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss); - GetNodalFunctions(basis,gauss); - - vx_input->GetInputValue(&vx,gauss); - - D_scalar=gauss->weight*Jdet; - - TripleMultiply(basis,1,numnodes,1, - &D_scalar,1,1,0, - basis,1,numnodes,0, - &Ke->values[0],1); - - GetNodalFunctions(B,gauss); - GetBprimeMasstransport(Bprime,&xyz_list[0][0],gauss); - - D_scalar=dt*gauss->weight*Jdet*vx; - TripleMultiply(B,1,numnodes,1, - &D_scalar,1,1,0, - Bprime,1,numnodes,0, - &Ke->values[0],1); - - if(stabilization==2){ - /*Streamline upwinding*/ - vel=fabs(vx)+1.e-8; - D_scalar=dt*gauss->weight*Jdet*h/(2.*vel)*vx; - } - else if(stabilization==1){ - /*SSA*/ - vx_input->GetInputAverage(&vx); - D_scalar=dt*gauss->weight*Jdet*h/2.*fabs(vx); - } - if(stabilization==1 || stabilization==2){ - TripleMultiply(Bprime,1,numnodes,1, - &D_scalar,1,1,0, - Bprime,1,numnodes,0, - &Ke->values[0],1); - } - } - - /*Clean up and return*/ - xDelete(basis); - xDelete(B); - xDelete(Bprime); - delete gauss; - return Ke; -} -/*}}}*/ -/*FUNCTION Seg::CreateMassMatrix {{{*/ -ElementMatrix* Seg::CreateMassMatrix(void){ - - /* Intermediaries */ - IssmDouble D,Jdet; - IssmDouble xyz_list[NUMVERTICES][3]; - - /*Fetch number of nodes and dof for this finite element*/ - int numnodes = this->NumberofNodes(); - - /*Initialize Element matrix and vectors*/ - ElementMatrix* Ke = new ElementMatrix(nodes,numnodes,this->parameters,NoneApproximationEnum); - IssmDouble* basis = xNew(numnodes); - - /*Retrieve all inputs and parameters*/ - ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); - - /* Start looping on the number of gaussian points: */ - GaussSeg* gauss=new GaussSeg(2); - for(int ig=gauss->begin();igend();ig++){ - - gauss->GaussPoint(ig); - - GetNodalFunctions(basis,gauss); - GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss); - D=gauss->weight*Jdet; - - TripleMultiply(basis,1,numnodes,1, - &D,1,1,0, - basis,1,numnodes,0, - &Ke->values[0],1); - } - - /*Clean up and return*/ - delete gauss; - xDelete(basis); - return Ke; -} -/*}}}*/ /*FUNCTION Seg::GetNumberOfNodes;{{{*/ int Seg::GetNumberOfNodes(void){ Index: /issm/trunk-jpl/src/c/classes/Elements/Seg.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Seg.h (revision 16992) +++ /issm/trunk-jpl/src/c/classes/Elements/Seg.h (revision 16993) @@ -69,7 +69,5 @@ void SetCurrentConfiguration(Elements* elements,Loads* loads,Nodes* nodes,Materials* materials,Parameters* parameters){_error_("not implemented yet");}; void SetwiseNodeConnectivity(int* d_nz,int* o_nz,Node* node,bool* flags,int* flagsindices,int set1_enum,int set2_enum){_error_("not implemented yet");}; - ElementMatrix* CreateKMatrix(void){_error_("not implemented yet");}; void CreateDVector(Vector* df){_error_("not implemented yet");}; - void CreateJacobianMatrix(Matrix* Jff){_error_("not implemented yet");}; void Delta18oParameterization(void){_error_("not implemented yet");}; void ElementSizes(IssmDouble* hx,IssmDouble* hy,IssmDouble* hz){_error_("not implemented yet");}; @@ -257,7 +255,4 @@ /*}}}*/ /*Seg specific routines:*/ - ElementMatrix* CreateMassMatrix(void); - ElementMatrix* CreateKMatrixFreeSurfaceTop(void); - ElementMatrix* CreateKMatrixFreeSurfaceBase(void); IssmDouble GetSize(void); }; Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 16992) +++ /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 16993) @@ -191,218 +191,8 @@ } /*}}}*/ -/*FUNCTION Tria::CreateKMatrix(void) {{{*/ -ElementMatrix* Tria::CreateKMatrix(void){ - - /*retreive parameters: */ - int analysis_type; - int meshtype; - parameters->FindParam(&analysis_type,AnalysisTypeEnum); - - /*Checks in debugging mode{{{*/ - _assert_(this->nodes && this->material && this->matpar && this->parameters && this->inputs); - /*}}}*/ - - /*Skip if water element*/ - if(NoIceInElement()) return NULL; - - /*Just branch to the correct element stiffness matrix generator, according to the type of analysis we are carrying out: */ - switch(analysis_type){ - #ifdef _HAVE_STRESSBALANCE_ - case StressbalanceAnalysisEnum: - int approximation; - inputs->GetInputValue(&approximation,ApproximationEnum); - switch(approximation){ - case SSAApproximationEnum: - return CreateKMatrixStressbalanceSSA(); - case FSApproximationEnum: - return CreateKMatrixStressbalanceFS(); - case SIAApproximationEnum: - return NULL; - case NoneApproximationEnum: - return NULL; - default: - _error_("Approximation " << EnumToStringx(approximation) << " not supported yet"); - } - break; - case StressbalanceSIAAnalysisEnum: - return CreateKMatrixStressbalanceSIA(); - break; - #endif - #ifdef _HAVE_DAMAGE_ - case DamageEvolutionAnalysisEnum: - return CreateKMatrixDamageEvolution(); - break; - #endif - case L2ProjectionBaseAnalysisEnum: - parameters->FindParam(&meshtype,MeshTypeEnum); - if(meshtype==Mesh2DverticalEnum){ - return CreateBasalMassMatrix(); - } - else{ - return CreateMassMatrix(); - } - break; - #ifdef _HAVE_MASSTRANSPORT_ - case MasstransportAnalysisEnum: - return CreateKMatrixMasstransport(); - break; - case FreeSurfaceTopAnalysisEnum: - parameters->FindParam(&meshtype,MeshTypeEnum); - switch(meshtype){ - case Mesh2DverticalEnum: - return CreateKMatrixFreeSurfaceTop1D(); - case Mesh3DEnum: - return CreateKMatrixFreeSurfaceTop(); - default: - _error_("Mesh not supported yet"); - } - break; - case FreeSurfaceBaseAnalysisEnum: - parameters->FindParam(&meshtype,MeshTypeEnum); - switch(meshtype){ - case Mesh2DverticalEnum: - return CreateKMatrixFreeSurfaceBase1D(); - case Mesh3DEnum: - return CreateKMatrixFreeSurfaceBase(); - default: - _error_("Mesh not supported yet"); - } - break; - case ExtrudeFromBaseAnalysisEnum: case ExtrudeFromTopAnalysisEnum: - return CreateKMatrixExtrusion(); - #endif - #ifdef _HAVE_HYDROLOGY_ - case HydrologyShreveAnalysisEnum: - return CreateKMatrixHydrologyShreve(); - break; - case HydrologyDCInefficientAnalysisEnum: - return CreateKMatrixHydrologyDCInefficient(); - break; - case HydrologyDCEfficientAnalysisEnum: - return CreateKMatrixHydrologyDCEfficient(); - break; - case L2ProjectionEPLAnalysisEnum: - return CreateEPLDomainMassMatrix(); - break; - #endif - #ifdef _HAVE_BALANCED_ - case BalancethicknessAnalysisEnum: - return CreateKMatrixBalancethickness(); - break; - case BalancevelocityAnalysisEnum: - return CreateKMatrixBalancevelocity(); - break; - case SmoothedSurfaceSlopeXAnalysisEnum: case SmoothedSurfaceSlopeYAnalysisEnum: - return CreateKMatrixSmoothedSlope(); - break; - #endif - #ifdef _HAVE_CONTROL_ - case AdjointBalancethicknessAnalysisEnum: - return CreateKMatrixAdjointBalancethickness(); - break; - case AdjointHorizAnalysisEnum: - return CreateKMatrixAdjointSSA(); - break; - #endif - default: - _error_("analysis " << analysis_type << " (" << EnumToStringx(analysis_type) << ") not supported yet"); - } - - /*Make compiler happy*/ - return NULL; -} -/*}}}*/ -/*FUNCTION Tria::CreateMassMatrix {{{*/ -ElementMatrix* Tria::CreateMassMatrix(void){ - - /* Intermediaries */ - IssmDouble D,Jdet; - IssmDouble xyz_list[NUMVERTICES][3]; - - /*Fetch number of nodes and dof for this finite element*/ - int numnodes = this->NumberofNodes(); - - /*Initialize Element matrix and vectors*/ - ElementMatrix* Ke = new ElementMatrix(nodes,numnodes,this->parameters,NoneApproximationEnum); - IssmDouble* basis = xNew(numnodes); - - /*Retrieve all inputs and parameters*/ - ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); - - /* Start looping on the number of gaussian points: */ - GaussTria* gauss=new GaussTria(2); - for(int ig=gauss->begin();igend();ig++){ - - gauss->GaussPoint(ig); - - GetNodalFunctions(basis,gauss); - GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss); - D=gauss->weight*Jdet; - - TripleMultiply(basis,1,numnodes,1, - &D,1,1,0, - basis,1,numnodes,0, - &Ke->values[0],1); - } - - /*Clean up and return*/ - delete gauss; - xDelete(basis); - return Ke; -} -/*}}}*/ -/*FUNCTION Tria::CreateBasalMassMatrix{{{*/ -ElementMatrix* Tria::CreateBasalMassMatrix(void){ - - if(!HasEdgeOnBed()) return NULL; - - int index1,index2; - this->EdgeOnBedIndices(&index1,&index2); - - Seg* seg=(Seg*)SpawnSeg(index1,index2); - ElementMatrix* Ke=seg->CreateMassMatrix(); - delete seg->material; delete seg; - - /*clean up and return*/ - return Ke; -} -/*}}}*/ /*FUNCTION Tria::CreateDVector {{{*/ void Tria::CreateDVector(Vector* df){ /*Nothing done yet*/ -} -/*}}}*/ -/*FUNCTION Tria::CreateJacobianMatrix{{{*/ -void Tria::CreateJacobianMatrix(Matrix* Jff){ - - /*retrieve parameters: */ - ElementMatrix* Ke=NULL; - int analysis_type; - parameters->FindParam(&analysis_type,AnalysisTypeEnum); - - /*Checks in debugging {{{*/ - _assert_(this->nodes && this->material && this->matpar && this->parameters && this->inputs); - /*}}}*/ - - /*Skip if water element*/ - if(NoIceInElement()) return; - - /*Just branch to the correct element stiffness matrix generator, according to the type of analysis we are carrying out: */ - switch(analysis_type){ -#ifdef _HAVE_STRESSBALANCE_ - case StressbalanceAnalysisEnum: - Ke=CreateJacobianStressbalanceSSA(); - break; -#endif - default: - _error_("analysis " << analysis_type << " (" << EnumToStringx(analysis_type) << ") not supported yet"); - } - - /*Add to global matrix*/ - if(Ke){ - Ke->AddToGlobal(Jff); - delete Ke; - } } /*}}}*/ @@ -3273,417 +3063,4 @@ #ifdef _HAVE_STRESSBALANCE_ -/*FUNCTION Tria::CreateKMatrixStressbalanceFS{{{*/ -ElementMatrix* Tria::CreateKMatrixStressbalanceFS(void){ - - ElementMatrix* Ke1 = NULL; - ElementMatrix* Ke2 = NULL; - ElementMatrix* Ke = NULL; - - /*compute all stiffness matrices for this element*/ - Ke1=CreateKMatrixStressbalanceFSViscous(); - Ke2=CreateKMatrixStressbalanceFSFriction(); - Ke =new ElementMatrix(Ke1,Ke2); - - /*clean-up and return*/ - delete Ke1; - delete Ke2; - return Ke; - -} -/*}}}*/ -/*FUNCTION Tria::CreateKMatrixStressbalanceFSViscous {{{*/ -ElementMatrix* Tria::CreateKMatrixStressbalanceFSViscous(void){ - - /*Intermediaries */ - int i,approximation; - IssmDouble Jdet,viscosity,FSreconditioning,D_scalar; - IssmDouble xyz_list[NUMVERTICES][3]; - IssmDouble epsilon[3]; /* epsilon=[exx,eyy,exy];*/ - GaussTria *gauss=NULL; - - /*Fetch number of nodes and dof for this finite element*/ - int vnumnodes = this->NumberofNodesVelocity(); - int pnumnodes = this->NumberofNodesPressure(); - int numdof = vnumnodes*NDOF2 + pnumnodes*NDOF1; - - /*Prepare coordinate system list*/ - int* cs_list = xNew(vnumnodes+pnumnodes); - for(i=0;iparameters,FSvelocityEnum); - IssmDouble* B = xNew(5*numdof); - IssmDouble* Bprime = xNew(5*numdof); - IssmDouble* D = xNewZeroInit(5*5); - - /*Retrieve all inputs and parameters*/ - ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); - parameters->FindParam(&FSreconditioning,StressbalanceFSreconditioningEnum); - Input* vx_input=inputs->GetInput(VxEnum); _assert_(vx_input); - Input* vy_input=inputs->GetInput(VyEnum); _assert_(vy_input); - - /* Start looping on the number of gaussian points: */ - gauss=new GaussTria(5); - for(int ig=gauss->begin();igend();ig++){ - - gauss->GaussPoint(ig); - - GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss); - GetBFS(B,&xyz_list[0][0],gauss); - GetBprimeFS(Bprime,&xyz_list[0][0],gauss); - - this->StrainRateSSA(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input); - material->GetViscosity2dvertical(&viscosity,&epsilon[0]); - - D_scalar=gauss->weight*Jdet; - for(i=0;i<3;i++) D[i*5+i] = +D_scalar*2.*viscosity; - for(i=3;i<5;i++) D[i*5+i] = -D_scalar*FSreconditioning; - - TripleMultiply(B,5,numdof,1, - D,5,5,0, - Bprime,5,numdof,0, - Ke->values,1); - } - - /*Transform Coordinate System*/ - ::TransformStiffnessMatrixCoord(Ke,nodes,(vnumnodes+pnumnodes),cs_list); - - /*Clean up and return*/ - delete gauss; - xDelete(cs_list); - xDelete(B); - xDelete(Bprime); - xDelete(D); - return Ke; -} -/*}}}*/ -/*FUNCTION Tria::CreateKMatrixStressbalanceFSFriction{{{*/ -ElementMatrix* Tria::CreateKMatrixStressbalanceFSFriction(void){ - - /*Intermediaries */ - int i,j; - int analysis_type; - int indices[2]; - IssmDouble alpha2,Jdet; - IssmDouble xyz_list[NUMVERTICES][3]; - IssmDouble xyz_list_seg[NUMVERTICES1D][3]; - Friction *friction = NULL; - - /*Initialize Element matrix and return if necessary*/ - if(IsFloating() || !HasEdgeOnBed()) return NULL; - - /*Fetch number of nodes and dof for this finite element*/ - int vnumnodes = this->NumberofNodesVelocity(); - int pnumnodes = this->NumberofNodesPressure(); - int numdof = vnumnodes*NDOF2 + pnumnodes*NDOF1; - - /*Prepare coordinate system list*/ - int* cs_list = xNew(vnumnodes+pnumnodes); - for(i=0;iparameters,FSvelocityEnum); - IssmDouble* BFriction = xNew(1*numdof); - IssmDouble D; - - /*Retrieve all inputs and parameters*/ - ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); - parameters->FindParam(&analysis_type,AnalysisTypeEnum); - Input* vx_input=inputs->GetInput(VxEnum); _assert_(vx_input); - Input* vy_input=inputs->GetInput(VyEnum); _assert_(vy_input); - - /*Get vertex indices that lie on bed*/ - this->EdgeOnBedIndices(&indices[0],&indices[1]); - for(i=0;ibegin();igend();ig++){ - - gauss->GaussPoint(ig); - - GetSegmentJacobianDeterminant(&Jdet,&xyz_list_seg[0][0],gauss); - GetLFS(BFriction,gauss); - - friction->GetAlpha2(&alpha2,gauss,vx_input,vy_input,NULL); - D = +alpha2*gauss->weight*Jdet; //taub_x = -alpha2 vx - - TripleMultiply(BFriction,1,numdof,1, - &D,1,1,0, - BFriction,1,numdof,0, - Ke->values,1); - } - - /*DO NOT Transform Coordinate System: this stiffness matrix is already expressed in tangential coordinates*/ - //::TransformStiffnessMatrixCoord(Ke,nodes,(vnumnodes+pnumnodes),cs_list); - - /*Clean up and return*/ - delete gauss; - delete friction; - xDelete(BFriction); - xDelete(cs_list); - return Ke; -} -/*}}}*/ -/*FUNCTION Tria::CreateKMatrixStressbalanceSSA {{{*/ -ElementMatrix* Tria::CreateKMatrixStressbalanceSSA(void){ - - /*compute all stiffness matrices for this element*/ - ElementMatrix* Ke1=CreateKMatrixStressbalanceSSAViscous(); - ElementMatrix* Ke2=CreateKMatrixStressbalanceSSAFriction(); - ElementMatrix* Ke =new ElementMatrix(Ke1,Ke2); - - /*clean-up and return*/ - delete Ke1; - delete Ke2; - return Ke; -} -/*}}}*/ -/*FUNCTION Tria::CreateKMatrixStressbalanceSSAViscous{{{*/ -ElementMatrix* Tria::CreateKMatrixStressbalanceSSAViscous(void){ - - /*Intermediaries*/ - IssmDouble xyz_list[NUMVERTICES][3]; - IssmDouble viscosity,newviscosity,oldviscosity; - IssmDouble viscosity_overshoot,thickness,Jdet; - IssmDouble epsilon[3],oldepsilon[3]; /* epsilon=[exx,eyy,exy]; */ - IssmDouble D_scalar; - - /*Fetch number of nodes and dof for this finite element*/ - int numnodes = this->NumberofNodes(); - int numdof = numnodes*NDOF2; - - /*Initialize Element matrix and vectors*/ - ElementMatrix* Ke = new ElementMatrix(nodes,numnodes,this->parameters,SSAApproximationEnum); - IssmDouble* B = xNew(3*numdof); - IssmDouble* Bprime = xNew(3*numdof); - IssmDouble* D = xNewZeroInit(3*3); - - /*Retrieve all inputs and parameters*/ - ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); - Input* thickness_input=inputs->GetInput(ThicknessEnum); _assert_(thickness_input); - Input* vx_input=inputs->GetInput(VxEnum); _assert_(vx_input); - Input* vy_input=inputs->GetInput(VyEnum); _assert_(vy_input); - Input* vxold_input=inputs->GetInput(VxPicardEnum); _assert_(vxold_input); - Input* vyold_input=inputs->GetInput(VyPicardEnum); _assert_(vyold_input); - this->parameters->FindParam(&viscosity_overshoot,StressbalanceViscosityOvershootEnum); - - /* Start looping on the number of gaussian points: */ - GaussTria* gauss = new GaussTria(2); - for(int ig=gauss->begin();igend();ig++){ - - gauss->GaussPoint(ig); - - GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss); - GetBSSA(&B[0], &xyz_list[0][0], gauss); - GetBprimeSSA(&Bprime[0], &xyz_list[0][0], gauss); - - this->StrainRateSSA(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input); - this->StrainRateSSA(&oldepsilon[0],&xyz_list[0][0],gauss,vxold_input,vyold_input); - material->GetViscosity2d(&viscosity, &epsilon[0]); - material->GetViscosity2d(&oldviscosity, &oldepsilon[0]); - thickness_input->GetInputValue(&thickness, gauss); - - newviscosity=viscosity+viscosity_overshoot*(viscosity-oldviscosity); - D_scalar=2.*newviscosity*thickness*gauss->weight*Jdet; - for(int i=0;i<3;i++) D[i*3+i]=D_scalar; - - TripleMultiply(B,3,numdof,1, - D,3,3,0, - Bprime,3,numdof,0, - &Ke->values[0],1); - } - - /*Transform Coordinate System*/ - ::TransformStiffnessMatrixCoord(Ke,nodes,numnodes,XYEnum); - - /*Clean up and return*/ - delete gauss; - xDelete(D); - xDelete(Bprime); - xDelete(B); - return Ke; -} -/*}}}*/ -/*FUNCTION Tria::CreateKMatrixStressbalanceSSAFriction {{{*/ -ElementMatrix* Tria::CreateKMatrixStressbalanceSSAFriction(void){ - - /*Intermediaries*/ - bool mainlyfloating; - int analysis_type,migration_style; - int point1; - IssmDouble alpha2,Jdet; - IssmDouble fraction1,fraction2; - IssmDouble phi=1.0; - IssmDouble D_scalar; - IssmDouble gllevelset; - IssmDouble xyz_list[NUMVERTICES][3]; - Friction *friction = NULL; - GaussTria *gauss = NULL; - - /*Return if element is inactive*/ - if(IsFloating()) return NULL; - - /*Fetch number of nodes and dof for this finite element*/ - int numnodes = this->NumberofNodes(); - int numdof = numnodes*NDOF2; - - /*Initialize Element matrix and vectors*/ - ElementMatrix* Ke = new ElementMatrix(nodes,numnodes,this->parameters,SSAApproximationEnum); - IssmDouble* B = xNew(2*numdof); - IssmDouble* D = xNewZeroInit(2*2); - - /*Retrieve all inputs and parameters*/ - ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); - parameters->FindParam(&migration_style,GroundinglineMigrationEnum); - Input* surface_input=inputs->GetInput(SurfaceEnum); _assert_(surface_input); - Input* vx_input=inputs->GetInput(VxEnum); _assert_(vx_input); - Input* vy_input=inputs->GetInput(VyEnum); _assert_(vy_input); - Input* gllevelset_input=NULL; - parameters->FindParam(&analysis_type,AnalysisTypeEnum); - - /*build friction object, used later on: */ - friction=new Friction(this,2); - - /*Recover portion of element that is grounded*/ - if(migration_style==SubelementMigrationEnum) phi=this->GetGroundedPortion(&xyz_list[0][0]); - if(migration_style==SubelementMigration2Enum){ - gllevelset_input=inputs->GetInput(MaskGroundediceLevelsetEnum); _assert_(gllevelset_input); - this->GetGroundedPart(&point1,&fraction1,&fraction2,&mainlyfloating); - gauss=new GaussTria(point1,fraction1,fraction2,mainlyfloating,2); - } - else{ - gauss=new GaussTria(2); - } - - /* Start looping on the number of gaussian points: */ - for(int ig=gauss->begin();igend();ig++){ - - gauss->GaussPoint(ig); - - friction->GetAlpha2(&alpha2, gauss,vx_input,vy_input,NULL); - if(migration_style==SubelementMigrationEnum) alpha2=phi*alpha2; - if(migration_style==SubelementMigration2Enum){ - gllevelset_input->GetInputValue(&gllevelset, gauss); - if(gllevelset<0) alpha2=0; - } - - GetBSSAFriction(&B[0], &xyz_list[0][0], gauss); - GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss); - D_scalar=alpha2*gauss->weight*Jdet; - for(int i=0;i<2;i++) D[i*2+i]=D_scalar; - - TripleMultiply(B,2,numdof,1, - D,2,2,0, - B,2,numdof,0, - &Ke->values[0],1); - } - - /*Transform Coordinate System*/ - ::TransformStiffnessMatrixCoord(Ke,nodes,numnodes,XYEnum); - - /*Clean up and return*/ - delete gauss; - delete friction; - xDelete(D); - xDelete(B); - return Ke; -} -/*}}}*/ -/*FUNCTION Tria::CreateKMatrixStressbalanceSIA{{{*/ -ElementMatrix* Tria::CreateKMatrixStressbalanceSIA(void){ - - /*Intermediaries*/ - IssmDouble connectivity; - - /*Fetch number of nodes and dof for this finite element*/ - int numnodes = this->NumberofNodes(); _assert_(numnodes==3); - int numdof = numnodes*NDOF2; - - /*Initialize Element matrix*/ - ElementMatrix* Ke=new ElementMatrix(nodes,numnodes,this->parameters,NoneApproximationEnum); - - /*Create Element matrix*/ - for(int i=0;iConnectivity(); - Ke->values[(2*i+0)*numdof+(2*i+0)]=1./connectivity; - Ke->values[(2*i+1)*numdof+(2*i+1)]=1./connectivity; - } - - /*Clean up and return*/ - return Ke; -} -/*}}}*/ -/*FUNCTION Tria::CreateJacobianStressbalanceSSA{{{*/ -ElementMatrix* Tria::CreateJacobianStressbalanceSSA(void){ - - /*Intermediaries */ - int i,j; - IssmDouble xyz_list[NUMVERTICES][3]; - IssmDouble Jdet,thickness; - IssmDouble eps1dotdphii,eps1dotdphij; - IssmDouble eps2dotdphii,eps2dotdphij; - IssmDouble mu_prime; - IssmDouble epsilon[3];/* epsilon=[exx,eyy,exy];*/ - IssmDouble eps1[2],eps2[2]; - GaussTria* gauss = NULL; - - /*Fetch number of nodes and dof for this finite element*/ - int numnodes = this->NumberofNodes(); - - /*Initialize Element matrix, vectors and Gaussian points*/ - ElementMatrix* Ke=CreateKMatrixStressbalanceSSA(); //Initialize Jacobian with regular SSA (first part of the Gateau derivative) - IssmDouble* dbasis = xNew(2*numnodes); - - /*Retrieve all inputs and parameters*/ - ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); - Input* vx_input=inputs->GetInput(VxEnum); _assert_(vx_input); - Input* vy_input=inputs->GetInput(VyEnum); _assert_(vy_input); - Input* thickness_input=inputs->GetInput(ThicknessEnum); _assert_(thickness_input); - - /* Start looping on the number of gaussian points: */ - gauss=new GaussTria(2); - for(int ig=gauss->begin();igend();ig++){ - - gauss->GaussPoint(ig); - - GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss); - GetNodalFunctionsDerivatives(dbasis,&xyz_list[0][0],gauss); - - thickness_input->GetInputValue(&thickness, gauss); - this->StrainRateSSA(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input); - material->GetViscosity2dDerivativeEpsSquare(&mu_prime,&epsilon[0]); - eps1[0]=2*epsilon[0]+epsilon[1]; eps2[0]=epsilon[2]; - eps1[1]=epsilon[2]; eps2[1]=epsilon[0]+2*epsilon[1]; - - for(i=0;ivalues[2*numnodes*(2*i+0)+2*j+0]+=gauss->weight*Jdet*2*mu_prime*thickness*eps1dotdphij*eps1dotdphii; - Ke->values[2*numnodes*(2*i+0)+2*j+1]+=gauss->weight*Jdet*2*mu_prime*thickness*eps2dotdphij*eps1dotdphii; - Ke->values[2*numnodes*(2*i+1)+2*j+0]+=gauss->weight*Jdet*2*mu_prime*thickness*eps1dotdphij*eps2dotdphii; - Ke->values[2*numnodes*(2*i+1)+2*j+1]+=gauss->weight*Jdet*2*mu_prime*thickness*eps2dotdphij*eps2dotdphii; - } - } - } - - /*Transform Coordinate System*/ - ::TransformStiffnessMatrixCoord(Ke,nodes,numnodes,XYEnum); - - /*Clean up and return*/ - xDelete(dbasis); - delete gauss; - return Ke; -} -/*}}}*/ /*FUNCTION Tria::GetYcoord {{{*/ IssmDouble Tria::GetYcoord(GaussTria* gauss){ @@ -5037,98 +4414,4 @@ } /*}}}*/ -/*FUNCTION Tria::CreateKMatrixAdjointBalancethickness {{{*/ -ElementMatrix* Tria::CreateKMatrixAdjointBalancethickness(void){ - - ElementMatrix* Ke=NULL; - - /*Get Element Matrix of the forward model*/ - switch(GetElementType()){ - case P1Enum: - Ke=CreateKMatrixBalancethickness_CG(); - break; - case P1DGEnum: - Ke=CreateKMatrixBalancethickness_DG(); - break; - default: - _error_("Element type " << EnumToStringx(GetElementType()) << " not supported yet"); - } - - /*Transpose and return Ke*/ - Ke->Transpose(); - return Ke; -} -/*}}}*/ -/*FUNCTION Tria::CreateKMatrixAdjointSSA{{{*/ -ElementMatrix* Tria::CreateKMatrixAdjointSSA(void){ - - /*Intermediaries */ - int i,j; - bool incomplete_adjoint; - IssmDouble xyz_list[NUMVERTICES][3]; - IssmDouble Jdet,thickness; - IssmDouble eps1dotdphii,eps1dotdphij; - IssmDouble eps2dotdphii,eps2dotdphij; - IssmDouble mu_prime; - IssmDouble epsilon[3];/* epsilon=[exx,eyy,exy];*/ - IssmDouble eps1[2],eps2[2]; - GaussTria *gauss=NULL; - - /*Fetch number of nodes and dof for this finite element*/ - int numnodes = this->NumberofNodes(); - - /*Initialize Jacobian with regular SSA (first part of the Gateau derivative)*/ - parameters->FindParam(&incomplete_adjoint,InversionIncompleteAdjointEnum); - ElementMatrix* Ke=CreateKMatrixStressbalanceSSA(); - if(incomplete_adjoint) return Ke; - - /*Retrieve all inputs and parameters*/ - ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); - Input* vx_input=inputs->GetInput(VxEnum); _assert_(vx_input); - Input* vy_input=inputs->GetInput(VyEnum); _assert_(vy_input); - Input* thickness_input=inputs->GetInput(ThicknessEnum); _assert_(thickness_input); - - /*Allocate dbasis*/ - IssmDouble* dbasis = xNew(2*numnodes); - - /* Start looping on the number of gaussian points: */ - gauss=new GaussTria(2); - for(int ig=gauss->begin();igend();ig++){ - - gauss->GaussPoint(ig); - - GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss); - GetNodalFunctionsDerivatives(dbasis,&xyz_list[0][0],gauss); - - thickness_input->GetInputValue(&thickness, gauss); - this->StrainRateSSA(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input); - material->GetViscosity2dDerivativeEpsSquare(&mu_prime,&epsilon[0]); - eps1[0]=2*epsilon[0]+epsilon[1]; eps2[0]=epsilon[2]; - eps1[1]=epsilon[2]; eps2[1]=epsilon[0]+2*epsilon[1]; - - for(i=0;ivalues[2*numnodes*(2*i+0)+2*j+0]+=gauss->weight*Jdet*2*mu_prime*thickness*eps1dotdphij*eps1dotdphii; - Ke->values[2*numnodes*(2*i+0)+2*j+1]+=gauss->weight*Jdet*2*mu_prime*thickness*eps2dotdphij*eps1dotdphii; - Ke->values[2*numnodes*(2*i+1)+2*j+0]+=gauss->weight*Jdet*2*mu_prime*thickness*eps1dotdphij*eps2dotdphii; - Ke->values[2*numnodes*(2*i+1)+2*j+1]+=gauss->weight*Jdet*2*mu_prime*thickness*eps2dotdphij*eps2dotdphii; - } - } - } - - /*Transform Coordinate System*/ - ::TransformStiffnessMatrixCoord(Ke,nodes,numnodes,XYEnum); - - /*Clean up and return*/ - delete gauss; - xDelete(dbasis); - //Ke->Transpose(); - return Ke; -} -/*}}}*/ /*FUNCTION Tria::GetVectorFromControlInputs{{{*/ void Tria::GetVectorFromControlInputs(Vector* vector,int control_enum,int control_index,const char* data){ @@ -5191,50 +4474,4 @@ ((ControlInput*)input)->SetInput(new_input); -} -/*}}}*/ -#endif - -#ifdef _HAVE_THERMAL_ -/*FUNCTION Tria::CreateKMatrixMelting {{{*/ -ElementMatrix* Tria::CreateKMatrixMelting(void){ - - /*Intermediaries */ - IssmDouble heatcapacity,latentheat; - IssmDouble Jdet,D_scalar; - IssmDouble xyz_list[NUMVERTICES][3]; - - /*Fetch number of nodes and dof for this finite element*/ - int numnodes = this->NumberofNodes(); - - /*Initialize Element vector and vectors*/ - ElementMatrix* Ke=new ElementMatrix(nodes,numnodes,this->parameters,NoneApproximationEnum); - IssmDouble* basis = xNew(numnodes); - - /*Retrieve all inputs and parameters*/ - ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); - latentheat = matpar->GetLatentHeat(); - heatcapacity = matpar->GetHeatCapacity(); - - /* Start looping on the number of gauss (nodes on the bedrock) */ - GaussTria* gauss=new GaussTria(2); - for(int ig=gauss->begin();igend();ig++){ - - gauss->GaussPoint(ig); - - GetNodalFunctions(&basis[0],gauss); - GetJacobianDeterminant(&Jdet, &xyz_list[0][0], gauss); - - D_scalar=latentheat/heatcapacity*gauss->weight*Jdet; - - TripleMultiply(basis,numnodes,1,0, - &D_scalar,1,1,0, - basis,1,numnodes,0, - &Ke->values[0],1); - } - - /*Clean up and return*/ - xDelete(basis); - delete gauss; - return Ke; } /*}}}*/ @@ -5294,234 +4531,4 @@ this->inputs->AddInput(new TriaInput(HydrologyWaterVxEnum,vx,P1Enum)); this->inputs->AddInput(new TriaInput(HydrologyWaterVyEnum,vy,P1Enum)); -} -/*}}}*/ -/*FUNCTION Tria::CreateKMatrixHydrologyShreve{{{*/ -ElementMatrix* Tria::CreateKMatrixHydrologyShreve(void){ - - /*Intermediaries */ - IssmDouble diffusivity; - IssmDouble Jdet,D_scalar,dt,h; - IssmDouble vx,vy,vel,dvxdx,dvydy; - IssmDouble dvx[2],dvy[2]; - IssmDouble xyz_list[NUMVERTICES][3]; - - /*Skip if water or ice shelf element*/ - if(NoIceInElement() || IsFloating()) return NULL; - - /*Fetch number of nodes and dof for this finite element*/ - int numnodes = this->NumberofNodes(); - - /*Initialize Element matrix and vectors*/ - ElementMatrix* Ke = new ElementMatrix(nodes,numnodes,this->parameters,NoneApproximationEnum); - IssmDouble* basis = xNew(numnodes); - IssmDouble* B = xNew(2*numnodes); - IssmDouble* Bprime = xNew(2*numnodes); - IssmDouble D[2][2]; - - /*Create water velocity vx and vy from current inputs*/ - CreateHydrologyWaterVelocityInput(); - - /*Retrieve all inputs and parameters*/ - ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); - this->parameters->FindParam(&dt,TimesteppingTimeStepEnum); - this->parameters->FindParam(&diffusivity,HydrologyshreveStabilizationEnum); - Input* vx_input=inputs->GetInput(HydrologyWaterVxEnum); _assert_(vx_input); - Input* vy_input=inputs->GetInput(HydrologyWaterVyEnum); _assert_(vy_input); - h=sqrt(2*this->GetArea()); - - /* Start looping on the number of gaussian points: */ - GaussTria* gauss=new GaussTria(2); - for(int ig=gauss->begin();igend();ig++){ - - gauss->GaussPoint(ig); - - GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss); - GetNodalFunctions(basis,gauss); - - vx_input->GetInputValue(&vx,gauss); - vy_input->GetInputValue(&vy,gauss); - vx_input->GetInputDerivativeValue(&dvx[0],&xyz_list[0][0],gauss); - vy_input->GetInputDerivativeValue(&dvy[0],&xyz_list[0][0],gauss); - - D_scalar=gauss->weight*Jdet; - - TripleMultiply(basis,1,numnodes,1, - &D_scalar,1,1,0, - basis,1,numnodes,0, - Ke->values,1); - - GetBMasstransport(B,&xyz_list[0][0], gauss); - GetBprimeMasstransport(Bprime,&xyz_list[0][0], gauss); - - dvxdx=dvx[0]; - dvydy=dvy[1]; - D_scalar=dt*gauss->weight*Jdet; - - D[0][0]=D_scalar*dvxdx; - D[0][1]=0.; - D[1][0]=0.; - D[1][1]=D_scalar*dvydy; - TripleMultiply(B,2,numnodes,1, - &D[0][0],2,2,0, - B,2,numnodes,0, - &Ke->values[0],1); - - D[0][0]=D_scalar*vx; - D[1][1]=D_scalar*vy; - TripleMultiply(B,2,numnodes,1, - &D[0][0],2,2,0, - Bprime,2,numnodes,0, - &Ke->values[0],1); - - /*Artificial diffusivity*/ - vel=sqrt(vx*vx+vy*vy); - D[0][0]=D_scalar*diffusivity*h/(2*vel)*vx*vx; - D[1][0]=D_scalar*diffusivity*h/(2*vel)*vy*vx; - D[0][1]=D_scalar*diffusivity*h/(2*vel)*vx*vy; - D[1][1]=D_scalar*diffusivity*h/(2*vel)*vy*vy; - TripleMultiply(Bprime,2,numnodes,1, - &D[0][0],2,2,0, - Bprime,2,numnodes,0, - &Ke->values[0],1); - } - - /*Clean up and return*/ - xDelete(basis); - xDelete(B); - xDelete(Bprime); - delete gauss; - return Ke; -} -/*}}}*/ -/*FUNCTION Tria::CreateKMatrixHydrologyDCInefficient{{{*/ -ElementMatrix* Tria::CreateKMatrixHydrologyDCInefficient(void){ - - /* Intermediaries */ - IssmDouble D_scalar,Jdet; - IssmDouble sediment_transmitivity,dt; - IssmDouble sediment_storing; - IssmDouble xyz_list[NUMVERTICES][3]; - - /*Fetch number of nodes and dof for this finite element*/ - int numnodes = this->NumberofNodes(); - - /*Initialize Element matrix and vectors*/ - ElementMatrix* Ke = new ElementMatrix(nodes,numnodes,this->parameters); - IssmDouble* B = xNew(5*numnodes); - IssmDouble* basis = xNew(numnodes); - IssmDouble D[2][2]; - - /*Retrieve all inputs and parameters*/ - ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); - this->parameters->FindParam(&dt,TimesteppingTimeStepEnum); - sediment_storing = matpar->GetSedimentStoring(); - sediment_transmitivity = matpar->GetSedimentTransmitivity(); - - /* Start looping on the number of gaussian points: */ - GaussTria* gauss=new GaussTria(2); - for(int ig=gauss->begin();igend();ig++){ - - gauss->GaussPoint(ig); - - GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss); - - /*Diffusivity*/ - D_scalar=sediment_transmitivity*gauss->weight*Jdet; - if(reCast(dt)) D_scalar=D_scalar*dt; - D[0][0]=D_scalar; D[0][1]=0.; - D[1][0]=0.; D[1][1]=D_scalar; - GetBHydro(B,&xyz_list[0][0],gauss); - TripleMultiply(B,2,numnodes,1, - &D[0][0],2,2,0, - B,2,numnodes,0, - &Ke->values[0],1); - - /*Transient*/ - if(reCast(dt)){ - GetNodalFunctions(&basis[0],gauss); - D_scalar=sediment_storing*gauss->weight*Jdet; - - TripleMultiply(basis,numnodes,1,0, - &D_scalar,1,1,0, - basis,1,numnodes,0, - &Ke->values[0],1); - } - } - /*Clean up and return*/ - xDelete(basis); - xDelete(B); - delete gauss; - return Ke; -} -/*}}}*/ -/*FUNCTION Tria::CreateKMatrixHydrologyDCEfficient{{{*/ -ElementMatrix* Tria::CreateKMatrixHydrologyDCEfficient(void){ - - /* Intermediaries */ - IssmDouble D_scalar,Jdet,dt; - IssmDouble epl_thickness; - IssmDouble epl_conductivity; - IssmDouble epl_specificstoring; - IssmDouble xyz_list[NUMVERTICES][3]; - - /*Check that all nodes are active, else return empty matrix*/ - if(!this->AllActive()){ - return NULL; - } - - /*Fetch number of nodes and dof for this finite element*/ - int numnodes = this->NumberofNodes(); - - /*Initialize Element matrix and vectors*/ - ElementMatrix* Ke = new ElementMatrix(nodes,numnodes,this->parameters); - IssmDouble* B = xNew(5*numnodes); - IssmDouble* basis = xNew(numnodes); - IssmDouble D[2][2]; - - /*Retrieve all inputs and parameters*/ - ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); - Input* thickness_input=inputs->GetInput(HydrologydcEplThicknessEnum); _assert_(thickness_input); - this->parameters->FindParam(&dt,TimesteppingTimeStepEnum); - epl_specificstoring = matpar->GetEplSpecificStoring(); - epl_conductivity = matpar->GetEplConductivity(); - - - /* Start looping on the number of gaussian points: */ - GaussTria* gauss=new GaussTria(2); - for(int ig=gauss->begin();igend();ig++){ - - - gauss->GaussPoint(ig); - GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss); - thickness_input->GetInputValue(&epl_thickness,gauss); - - /*Diffusivity*/ - D_scalar=epl_conductivity*epl_thickness*gauss->weight*Jdet; - if(reCast(dt)) D_scalar=D_scalar*dt; - D[0][0]=D_scalar; D[0][1]=0.; - D[1][0]=0.; D[1][1]=D_scalar; - GetBHydro(B,&xyz_list[0][0],gauss); - TripleMultiply(B,2,numnodes,1, - &D[0][0],2,2,0, - B,2,numnodes,0, - &Ke->values[0],1); - - /*Transient*/ - if(reCast(dt)){ - GetNodalFunctions(basis,gauss); - D_scalar=epl_specificstoring*epl_thickness*gauss->weight*Jdet; - - TripleMultiply(basis,numnodes,1,0, - &D_scalar,1,1,0, - basis,1,numnodes,0, - &Ke->values[0],1); - } - } - - /*Clean up and return*/ - xDelete(basis); - xDelete(B); - delete gauss; - return Ke; } /*}}}*/ @@ -5898,730 +4905,5 @@ #endif -#ifdef _HAVE_MASSTRANSPORT_ -/*FUNCTION Tria::CreateKMatrixExtrusion {{{*/ -ElementMatrix* Tria::CreateKMatrixExtrusion(void){ - - /*compute all stiffness matrices for this element*/ - ElementMatrix* Ke1=CreateKMatrixExtrusionVolume(); - ElementMatrix* Ke2=CreateKMatrixExtrusionSurface(); - ElementMatrix* Ke3=CreateKMatrixExtrusionBed(); - ElementMatrix* Ke =new ElementMatrix(Ke1,Ke2,Ke3); - - /*clean-up and return*/ - delete Ke1; - delete Ke2; - delete Ke3; - return Ke; - -} -/*}}}*/ -/*FUNCTION Tria::CreateKMatrixExtrusionVolume {{{*/ -ElementMatrix* Tria::CreateKMatrixExtrusionVolume(void){ - - /*Intermediaries */ - IssmDouble Jdet; - IssmDouble xyz_list[NUMVERTICES][3]; - IssmDouble B[NDOF1][NUMVERTICES]; - IssmDouble Bprime[NDOF1][NUMVERTICES]; - IssmDouble DL_scalar; - GaussTria *gauss=NULL; - - /*Initialize Element matrix*/ - ElementMatrix* Ke=new ElementMatrix(nodes,NUMVERTICES,this->parameters,NoneApproximationEnum); - - /*Retrieve all inputs and parameters*/ - ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); - - /* Start looping on the number of gaussian points: */ - gauss=new GaussTria(2); - for(int ig=gauss->begin();igend();ig++){ - - gauss->GaussPoint(ig); - - GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss); - GetBExtrusion(&B[0][0], &xyz_list[0][0], gauss); - GetNodalFunctions(&Bprime[0][0],gauss); - - DL_scalar=gauss->weight*Jdet; - - TripleMultiply(&B[0][0],1,NUMVERTICES,1, - &DL_scalar,1,1,0, - &Bprime[0][0],1,NUMVERTICES,0, - &Ke->values[0],1); - } - - /*Clean up and return*/ - delete gauss; - return Ke; -} -/*}}}*/ -/*FUNCTION Tria::CreateKMatrixExtrusionSurface {{{*/ -ElementMatrix* Tria::CreateKMatrixExtrusionSurface(void){ - - if (!HasEdgeOnSurface()) return NULL; - - /*Constants*/ - const int numdof=NDOF1*NUMVERTICES; - - /*Intermediaries */ - int indices[2]; - IssmDouble xyz_list[NUMVERTICES][3]; - IssmDouble xyz_list_seg[NUMVERTICES1D][3]; - IssmDouble normal[3]; - IssmDouble Jdet,DL_scalar; - IssmDouble basis[NUMVERTICES]; - GaussTria *gauss=NULL; - - /*Initialize Element matrix*/ - ElementMatrix* Ke=new ElementMatrix(nodes,NUMVERTICES,this->parameters,NoneApproximationEnum); - - /*Retrieve all inputs and parameters*/ - ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); - - /*Get vertex indices that lie on bed*/ - this->EdgeOnSurfaceIndices(&indices[0],&indices[1]); - for(int i=0;ibegin();igend();ig++){ - - gauss->GaussPoint(ig); - - GetSegmentJacobianDeterminant(&Jdet,&xyz_list_seg[0][0],gauss); - GetNodalFunctions(&basis[0], gauss); - - DL_scalar= - gauss->weight*Jdet*normal[1]; - - TripleMultiply( basis,1,numdof,1, - &DL_scalar,1,1,0, - basis,1,numdof,0, - &Ke->values[0],1); - } - - /*Clean up and return*/ - delete gauss; - return Ke; -} -/*}}}*/ -/*FUNCTION Tria::CreateKMatrixExtrusionBed {{{*/ -ElementMatrix* Tria::CreateKMatrixExtrusionBed(void){ - - if (!HasEdgeOnBed()) return NULL; - - /*Constants*/ - const int numdof=NDOF1*NUMVERTICES; - - /*Intermediaries */ - int indices[2]; - IssmDouble xyz_list[NUMVERTICES][3]; - IssmDouble xyz_list_seg[NUMVERTICES1D][3]; - IssmDouble normal[3]; - IssmDouble Jdet,DL_scalar; - IssmDouble basis[NUMVERTICES]; - GaussTria *gauss=NULL; - - /*Initialize Element matrix*/ - ElementMatrix* Ke=new ElementMatrix(nodes,NUMVERTICES,this->parameters,NoneApproximationEnum); - - /*Retrieve all inputs and parameters*/ - ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); - - /*Get vertex indices that lie on bed*/ - this->EdgeOnBedIndices(&indices[0],&indices[1]); - for(int i=0;ibegin();igend();ig++){ - - gauss->GaussPoint(ig); - - GetSegmentJacobianDeterminant(&Jdet,&xyz_list_seg[0][0],gauss); - GetNodalFunctions(&basis[0], gauss); - - DL_scalar= - gauss->weight*Jdet*normal[1]; - - TripleMultiply( basis,1,numdof,1, - &DL_scalar,1,1,0, - basis,1,numdof,0, - &Ke->values[0],1); - } - - /*Clean up and return*/ - delete gauss; - return Ke; -} -/*}}}*/ -/*FUNCTION Tria::CreateKMatrixMasstransport {{{*/ -ElementMatrix* Tria::CreateKMatrixMasstransport(void){ - - switch(GetElementType()){ - case P1Enum: case P2Enum: - return CreateKMatrixMasstransport_CG(); - case P1DGEnum: - return CreateKMatrixMasstransport_DG(); - default: - _error_("Element type " << EnumToStringx(GetElementType()) << " not supported yet"); - } -} -/*}}}*/ -/*FUNCTION Tria::CreateKMatrixMasstransport_CG {{{*/ -ElementMatrix* Tria::CreateKMatrixMasstransport_CG(void){ - - /*Intermediaries */ - int stabilization; - int meshtype; - IssmDouble Jdet,D_scalar,dt,h; - IssmDouble vel,vx,vy,dvxdx,dvydy; - IssmDouble dvx[2],dvy[2]; - IssmDouble xyz_list[NUMVERTICES][3]; - - /*Fetch number of nodes for this finite element*/ - int numnodes = this->NumberofNodes(); - - /*Initialize Element matrix and vectors*/ - ElementMatrix* Ke = new ElementMatrix(nodes,numnodes,this->parameters,NoneApproximationEnum); - IssmDouble* basis = xNew(numnodes); - IssmDouble* B = xNew(2*numnodes); - IssmDouble* Bprime = xNew(2*numnodes); - IssmDouble D[2][2]; - - /*Retrieve all inputs and parameters*/ - ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); - this->parameters->FindParam(&dt,TimesteppingTimeStepEnum); - this->parameters->FindParam(&meshtype,MeshTypeEnum); - this->parameters->FindParam(&stabilization,MasstransportStabilizationEnum); - Input* vxaverage_input=NULL; - Input* vyaverage_input=NULL; - if(meshtype==Mesh2DhorizontalEnum){ - vxaverage_input=inputs->GetInput(VxEnum); _assert_(vxaverage_input); - vyaverage_input=inputs->GetInput(VyEnum); _assert_(vyaverage_input); - } - else{ - vxaverage_input=inputs->GetInput(VxAverageEnum); _assert_(vxaverage_input); - vyaverage_input=inputs->GetInput(VyAverageEnum); _assert_(vyaverage_input); - } - h=sqrt(2*this->GetArea()); - - /* Start looping on the number of gaussian points: */ - GaussTria *gauss=new GaussTria(2); - for(int ig=gauss->begin();igend();ig++){ - - gauss->GaussPoint(ig); - - GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss); - GetNodalFunctions(basis,gauss); - - vxaverage_input->GetInputValue(&vx,gauss); - vyaverage_input->GetInputValue(&vy,gauss); - vxaverage_input->GetInputDerivativeValue(&dvx[0],&xyz_list[0][0],gauss); - vyaverage_input->GetInputDerivativeValue(&dvy[0],&xyz_list[0][0],gauss); - - D_scalar=gauss->weight*Jdet; - - TripleMultiply(basis,1,numnodes,1, - &D_scalar,1,1,0, - basis,1,numnodes,0, - &Ke->values[0],1); - - GetBMasstransport(B,&xyz_list[0][0],gauss); - GetBprimeMasstransport(Bprime,&xyz_list[0][0],gauss); - - dvxdx=dvx[0]; - dvydy=dvy[1]; - D_scalar=dt*gauss->weight*Jdet; - - D[0][0]=D_scalar*dvxdx; - D[0][1]=0.; - D[1][0]=0.; - D[1][1]=D_scalar*dvydy; - TripleMultiply(B,2,numnodes,1, - &D[0][0],2,2,0, - B,2,numnodes,0, - &Ke->values[0],1); - - D[0][0]=D_scalar*vx; - D[1][1]=D_scalar*vy; - TripleMultiply(B,2,numnodes,1, - &D[0][0],2,2,0, - Bprime,2,numnodes,0, - &Ke->values[0],1); - - if(stabilization==2){ - /*Streamline upwinding*/ - vel=sqrt(vx*vx+vy*vy)+1.e-8; - D[0][0]=h/(2*vel)*vx*vx; - D[1][0]=h/(2*vel)*vy*vx; - D[0][1]=h/(2*vel)*vx*vy; - D[1][1]=h/(2*vel)*vy*vy; - } - else if(stabilization==1){ - /*SSA*/ - vxaverage_input->GetInputAverage(&vx); - vyaverage_input->GetInputAverage(&vy); - D[0][0]=h/2.0*fabs(vx); - D[0][1]=0.; - D[1][0]=0.; - D[1][1]=h/2.0*fabs(vy); - } - if(stabilization==1 || stabilization==2){ - D[0][0]=D_scalar*D[0][0]; - D[1][0]=D_scalar*D[1][0]; - D[0][1]=D_scalar*D[0][1]; - D[1][1]=D_scalar*D[1][1]; - TripleMultiply(Bprime,2,numnodes,1, - &D[0][0],2,2,0, - Bprime,2,numnodes,0, - &Ke->values[0],1); - } - } - - /*Clean up and return*/ - xDelete(basis); - xDelete(B); - xDelete(Bprime); - delete gauss; - return Ke; -} -/*}}}*/ -/*FUNCTION Tria::CreateKMatrixMasstransport_DG {{{*/ -ElementMatrix* Tria::CreateKMatrixMasstransport_DG(void){ - - /*Intermediaries */ - int meshtype; - IssmDouble xyz_list[NUMVERTICES][3]; - IssmDouble Jdet,D_scalar,dt,vx,vy; - - /*Fetch number of nodes for this finite element*/ - int numnodes = this->NumberofNodes(); - - /*Initialize Element matrix and vectors*/ - ElementMatrix* Ke = new ElementMatrix(nodes,numnodes,this->parameters,NoneApproximationEnum); - IssmDouble* basis = xNew(numnodes); - IssmDouble* B = xNew(2*numnodes); - IssmDouble* Bprime = xNew(2*numnodes); - IssmDouble D[2][2]; - - /*Retrieve all inputs and parameters*/ - ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); - this->parameters->FindParam(&dt,TimesteppingTimeStepEnum); - this->parameters->FindParam(&meshtype,MeshTypeEnum); - Input* vxaverage_input=NULL; - Input* vyaverage_input=NULL; - if(meshtype==Mesh2DhorizontalEnum){ - vxaverage_input=inputs->GetInput(VxEnum); _assert_(vxaverage_input); - vyaverage_input=inputs->GetInput(VyEnum); _assert_(vyaverage_input); - } - else{ - vxaverage_input=inputs->GetInput(VxAverageEnum); _assert_(vxaverage_input); - vyaverage_input=inputs->GetInput(VyAverageEnum); _assert_(vyaverage_input); - } - - /* Start looping on the number of gaussian points: */ - GaussTria* gauss=new GaussTria(2); - for(int ig=gauss->begin();igend();ig++){ - - gauss->GaussPoint(ig); - - vxaverage_input->GetInputValue(&vx,gauss); - vyaverage_input->GetInputValue(&vy,gauss); - - GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss); - GetNodalFunctions(basis,gauss); - - D_scalar=gauss->weight*Jdet; - - TripleMultiply(basis,1,numnodes,1, - &D_scalar,1,1,0, - basis,1,numnodes,0, - &Ke->values[0],1); - - /*WARNING: B and Bprime are inverted compared to usual masstransport!!!!*/ - GetBMasstransport(Bprime, &xyz_list[0][0], gauss); - GetBprimeMasstransport(B, &xyz_list[0][0], gauss); - - D_scalar=-dt*gauss->weight*Jdet; - D[0][0]=D_scalar*vx; - D[0][1]=0.; - D[1][0]=0.; - D[1][1]=D_scalar*vy; - - TripleMultiply(B,2,numnodes,1, - &D[0][0],2,2,0, - Bprime,2,numnodes,0, - &Ke->values[0],1); - } - - /*Clean up and return*/ - xDelete(basis); - xDelete(B); - xDelete(Bprime); - delete gauss; - return Ke; -} -/*}}}*/ -/*FUNCTION Tria::CreateKMatrixFreeSurfaceTop {{{*/ -ElementMatrix* Tria::CreateKMatrixFreeSurfaceTop(void){ - - /*Intermediaries */ - int stabilization; - IssmDouble Jdet,D_scalar,dt,h; - IssmDouble vel,vx,vy; - IssmDouble xyz_list[NUMVERTICES][3]; - - /*Fetch number of nodes for this finite element*/ - int numnodes = this->NumberofNodes(); - - /*Initialize Element matrix and vectors*/ - ElementMatrix* Ke = new ElementMatrix(nodes,numnodes,this->parameters,NoneApproximationEnum); - IssmDouble* basis = xNew(numnodes); - IssmDouble* B = xNew(2*numnodes); - IssmDouble* Bprime = xNew(2*numnodes); - IssmDouble D[2][2]; - - /*Retrieve all inputs and parameters*/ - ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); - this->parameters->FindParam(&dt,TimesteppingTimeStepEnum); - this->parameters->FindParam(&stabilization,MasstransportStabilizationEnum); - Input* vx_input=inputs->GetInput(VxEnum); _assert_(vx_input); - Input* vy_input=inputs->GetInput(VyEnum); _assert_(vy_input); - h=sqrt(2*this->GetArea()); - - /* Start looping on the number of gaussian points: */ - GaussTria *gauss=new GaussTria(2); - for(int ig=gauss->begin();igend();ig++){ - - gauss->GaussPoint(ig); - - GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss); - GetNodalFunctions(basis,gauss); - - vx_input->GetInputValue(&vx,gauss); - vy_input->GetInputValue(&vy,gauss); - - D_scalar=gauss->weight*Jdet; - - TripleMultiply(basis,1,numnodes,1, - &D_scalar,1,1,0, - basis,1,numnodes,0, - &Ke->values[0],1); - - GetBMasstransport(B,&xyz_list[0][0],gauss); - GetBprimeMasstransport(Bprime,&xyz_list[0][0],gauss); - - D_scalar=dt*gauss->weight*Jdet; - - D[0][0]=D_scalar*vx; - D[0][1]=0.; - D[1][0]=0.; - D[1][1]=D_scalar*vy; - TripleMultiply(B,2,numnodes,1, - &D[0][0],2,2,0, - Bprime,2,numnodes,0, - &Ke->values[0],1); - - if(stabilization==2){ - /*Streamline upwinding*/ - vel=sqrt(vx*vx+vy*vy)+1.e-8; - D[0][0]=h/(2*vel)*vx*vx; - D[1][0]=h/(2*vel)*vy*vx; - D[0][1]=h/(2*vel)*vx*vy; - D[1][1]=h/(2*vel)*vy*vy; - } - else if(stabilization==1){ - /*SSA*/ - vx_input->GetInputAverage(&vx); - vy_input->GetInputAverage(&vy); - D[0][0]=h/2.0*fabs(vx); - D[0][1]=0.; - D[1][0]=0.; - D[1][1]=h/2.0*fabs(vy); - } - if(stabilization==1 || stabilization==2){ - D[0][0]=D_scalar*D[0][0]; - D[1][0]=D_scalar*D[1][0]; - D[0][1]=D_scalar*D[0][1]; - D[1][1]=D_scalar*D[1][1]; - TripleMultiply(Bprime,2,numnodes,1, - &D[0][0],2,2,0, - Bprime,2,numnodes,0, - &Ke->values[0],1); - } - } - - /*Clean up and return*/ - xDelete(basis); - xDelete(B); - xDelete(Bprime); - delete gauss; - return Ke; -} -/*}}}*/ -/*FUNCTION Tria::CreateKMatrixFreeSurfaceTop1D {{{*/ -ElementMatrix* Tria::CreateKMatrixFreeSurfaceTop1D(void){ - - if(!HasEdgeOnSurface()) return NULL; - - int index1,index2; - this->EdgeOnSurfaceIndices(&index1,&index2); - - Seg* seg=(Seg*)SpawnSeg(index1,index2); - ElementMatrix* Ke=seg->CreateKMatrixFreeSurfaceTop(); - delete seg->material; delete seg; - - /*clean up and return*/ - return Ke; -} -/*}}}*/ -/*FUNCTION Tria::CreateKMatrixFreeSurfaceBase {{{*/ -ElementMatrix* Tria::CreateKMatrixFreeSurfaceBase(void){ - - /*Intermediaries */ - int stabilization; - IssmDouble Jdet,D_scalar,dt,h; - IssmDouble vel,vx,vy; - IssmDouble xyz_list[NUMVERTICES][3]; - - /*Fetch number of nodes for this finite element*/ - int numnodes = this->NumberofNodes(); - - /*Initialize Element matrix and vectors*/ - ElementMatrix* Ke = new ElementMatrix(nodes,numnodes,this->parameters,NoneApproximationEnum); - IssmDouble* basis = xNew(numnodes); - IssmDouble* B = xNew(2*numnodes); - IssmDouble* Bprime = xNew(2*numnodes); - IssmDouble D[2][2]; - - /*Retrieve all inputs and parameters*/ - ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); - this->parameters->FindParam(&dt,TimesteppingTimeStepEnum); - this->parameters->FindParam(&stabilization,MasstransportStabilizationEnum); - Input* vx_input=inputs->GetInput(VxEnum); _assert_(vx_input); - Input* vy_input=inputs->GetInput(VyEnum); _assert_(vy_input); - h=sqrt(2*this->GetArea()); - - /* Start looping on the number of gaussian points: */ - GaussTria *gauss=new GaussTria(2); - for(int ig=gauss->begin();igend();ig++){ - - gauss->GaussPoint(ig); - - GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss); - GetNodalFunctions(basis,gauss); - - vx_input->GetInputValue(&vx,gauss); - vy_input->GetInputValue(&vy,gauss); - - D_scalar=gauss->weight*Jdet; - - TripleMultiply(basis,1,numnodes,1, - &D_scalar,1,1,0, - basis,1,numnodes,0, - &Ke->values[0],1); - - GetBMasstransport(B,&xyz_list[0][0],gauss); - GetBprimeMasstransport(Bprime,&xyz_list[0][0],gauss); - - D_scalar=dt*gauss->weight*Jdet; - - D[0][0]=D_scalar*vx; - D[0][1]=0.; - D[1][0]=0.; - D[1][1]=D_scalar*vy; - TripleMultiply(B,2,numnodes,1, - &D[0][0],2,2,0, - Bprime,2,numnodes,0, - &Ke->values[0],1); - - if(stabilization==2){ - /*Streamline upwinding*/ - vel=sqrt(vx*vx+vy*vy)+1.e-8; - D[0][0]=h/(2*vel)*vx*vx; - D[1][0]=h/(2*vel)*vy*vx; - D[0][1]=h/(2*vel)*vx*vy; - D[1][1]=h/(2*vel)*vy*vy; - } - else if(stabilization==1){ - /*SSA*/ - vx_input->GetInputAverage(&vx); - vy_input->GetInputAverage(&vy); - D[0][0]=h/2.0*fabs(vx); - D[0][1]=0.; - D[1][0]=0.; - D[1][1]=h/2.0*fabs(vy); - } - if(stabilization==1 || stabilization==2){ - D[0][0]=D_scalar*D[0][0]; - D[1][0]=D_scalar*D[1][0]; - D[0][1]=D_scalar*D[0][1]; - D[1][1]=D_scalar*D[1][1]; - TripleMultiply(Bprime,2,numnodes,1, - &D[0][0],2,2,0, - Bprime,2,numnodes,0, - &Ke->values[0],1); - } - } - - /*Clean up and return*/ - xDelete(basis); - xDelete(B); - xDelete(Bprime); - delete gauss; - return Ke; -} -/*}}}*/ -/*FUNCTION Tria::CreateKMatrixFreeSurfaceBase1D {{{*/ -ElementMatrix* Tria::CreateKMatrixFreeSurfaceBase1D(void){ - - if(!HasEdgeOnBed()) return NULL; - - int index1,index2; - this->EdgeOnBedIndices(&index1,&index2); - - Seg* seg=(Seg*)SpawnSeg(index1,index2); - ElementMatrix* Ke=seg->CreateKMatrixFreeSurfaceBase(); - delete seg->material; delete seg; - - /*clean up and return*/ - return Ke; -} -/*}}}*/ -#endif - #ifdef _HAVE_DAMAGE_ -/*FUNCTION Tria::CreateKMatrixDamageEvolution {{{*/ -ElementMatrix* Tria::CreateKMatrixDamageEvolution(void){ - - switch(GetElementType()){ - case P1Enum: case P2Enum: - return CreateKMatrixDamageEvolution_CG(); - case P1DGEnum: - _error_("DG not implemented yet!");break; - default: - _error_("Element type " << EnumToStringx(GetElementType()) << " not supported yet"); - } -} -/*}}}*/ -/*FUNCTION Tria::CreateKMatrixDamageEvolution_CG {{{*/ -ElementMatrix* Tria::CreateKMatrixDamageEvolution_CG(void){ - - /*Intermediaries */ - int stabilization; - int meshtype; - IssmDouble Jdet,D_scalar,dt,h; - IssmDouble vel,vx,vy,dvxdx,dvydy; - IssmDouble dvx[2],dvy[2]; - IssmDouble xyz_list[NUMVERTICES][3]; - - /*Fetch number of nodes for this finite element*/ - int numnodes = this->NumberofNodes(); - - /*Initialize Element matrix and vectors*/ - ElementMatrix* Ke = new ElementMatrix(nodes,numnodes,this->parameters,NoneApproximationEnum); - IssmDouble* basis = xNew(numnodes); - IssmDouble* B = xNew(2*numnodes); - IssmDouble* Bprime = xNew(2*numnodes); - IssmDouble D[2][2]; - - /*Retrieve all inputs and parameters*/ - ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); - this->parameters->FindParam(&dt,TimesteppingTimeStepEnum); - this->parameters->FindParam(&meshtype,MeshTypeEnum); - this->parameters->FindParam(&stabilization,DamageStabilizationEnum); - Input* vxaverage_input=NULL; - Input* vyaverage_input=NULL; - if(meshtype==Mesh2DhorizontalEnum){ - vxaverage_input=inputs->GetInput(VxEnum); _assert_(vxaverage_input); - vyaverage_input=inputs->GetInput(VyEnum); _assert_(vyaverage_input); - } - else{ - vxaverage_input=inputs->GetInput(VxAverageEnum); _assert_(vxaverage_input); - vyaverage_input=inputs->GetInput(VyAverageEnum); _assert_(vyaverage_input); - } - h=sqrt(2*this->GetArea()); - - /* Start looping on the number of gaussian points: */ - GaussTria *gauss=new GaussTria(2); - for(int ig=gauss->begin();igend();ig++){ - - gauss->GaussPoint(ig); - - GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss); - GetNodalFunctions(basis,gauss); - - vxaverage_input->GetInputValue(&vx,gauss); - vyaverage_input->GetInputValue(&vy,gauss); - vxaverage_input->GetInputDerivativeValue(&dvx[0],&xyz_list[0][0],gauss); - vyaverage_input->GetInputDerivativeValue(&dvy[0],&xyz_list[0][0],gauss); - - D_scalar=gauss->weight*Jdet; - - TripleMultiply(basis,1,numnodes,1, - &D_scalar,1,1,0, - basis,1,numnodes,0, - &Ke->values[0],1); - GetBMasstransport(B,&xyz_list[0][0],gauss); - GetBprimeMasstransport(Bprime,&xyz_list[0][0],gauss); - - dvxdx=dvx[0]; - dvydy=dvy[1]; - D_scalar=dt*gauss->weight*Jdet; - - D[0][0]=D_scalar*dvxdx; - D[0][1]=0.; - D[1][0]=0.; - D[1][1]=D_scalar*dvydy; - TripleMultiply(B,2,numnodes,1, - &D[0][0],2,2,0, - B,2,numnodes,0, - &Ke->values[0],1); - - D[0][0]=D_scalar*vx; - D[1][1]=D_scalar*vy; - TripleMultiply(B,2,numnodes,1, - &D[0][0],2,2,0, - Bprime,2,numnodes,0, - &Ke->values[0],1); - - if(stabilization==2){ - /*Streamline upwinding*/ - vel=sqrt(vx*vx+vy*vy)+1.e-8; - D[0][0]=h/(2*vel)*vx*vx; - D[1][0]=h/(2*vel)*vy*vx; - D[0][1]=h/(2*vel)*vx*vy; - D[1][1]=h/(2*vel)*vy*vy; - } - else if(stabilization==1){ - /*SSA*/ - vxaverage_input->GetInputAverage(&vx); - vyaverage_input->GetInputAverage(&vy); - D[0][0]=h/2.0*fabs(vx); - D[0][1]=0.; - D[1][0]=0.; - D[1][1]=h/2.0*fabs(vy); - } - if(stabilization==1 || stabilization==2){ - D[0][0]=D_scalar*D[0][0]; - D[1][0]=D_scalar*D[1][0]; - D[0][1]=D_scalar*D[0][1]; - D[1][1]=D_scalar*D[1][1]; - TripleMultiply(Bprime,2,numnodes,1, - &D[0][0],2,2,0, - Bprime,2,numnodes,0, - &Ke->values[0],1); - } - } - - /*Clean up and return*/ - xDelete(basis); - xDelete(B); - xDelete(Bprime); - delete gauss; - return Ke; -} -/*}}}*/ /*FUNCTION Tria::DamageEvolutionF{{{*/ void Tria::DamageEvolutionF(IssmDouble* f){ @@ -6819,315 +5101,4 @@ } -} -/*}}}*/ -#endif - -#ifdef _HAVE_BALANCED_ -/*FUNCTION Tria::CreateKMatrixBalancethickness {{{*/ -ElementMatrix* Tria::CreateKMatrixBalancethickness(void){ - - switch(GetElementType()){ - case P1Enum: - return CreateKMatrixBalancethickness_CG(); - case P1DGEnum: - return CreateKMatrixBalancethickness_DG(); - default: - _error_("Element type " << EnumToStringx(GetElementType()) << " not supported yet"); - } - -} -/*}}}*/ -/*FUNCTION Tria::CreateKMatrixBalancethickness_CG {{{*/ -ElementMatrix* Tria::CreateKMatrixBalancethickness_CG(void){ - - /*Intermediaries */ - int stabilization,meshtype; - IssmDouble Jdet,vx,vy,dvxdx,dvydy,vel,h; - IssmDouble D_scalar; - IssmDouble dvx[2],dvy[2]; - IssmDouble xyz_list[NUMVERTICES][3]; - - /*Fetch number of nodes and dof for this finite element*/ - int numnodes = this->NumberofNodes(); - - /*Initialize Element matrix and vectors*/ - ElementMatrix* Ke = new ElementMatrix(nodes,numnodes,this->parameters,NoneApproximationEnum); - IssmDouble* basis = xNew(numnodes); - IssmDouble* B = xNew(2*numnodes); - IssmDouble* Bprime = xNew(2*numnodes); - IssmDouble D[2][2]; - - /*Retrieve all Inputs and parameters: */ - ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); - this->parameters->FindParam(&stabilization,BalancethicknessStabilizationEnum); - this->parameters->FindParam(&meshtype,MeshTypeEnum); - Input* vxaverage_input=NULL; - Input* vyaverage_input=NULL; - if(meshtype==Mesh2DhorizontalEnum){ - vxaverage_input=inputs->GetInput(VxEnum); _assert_(vxaverage_input); - vyaverage_input=inputs->GetInput(VyEnum); _assert_(vyaverage_input); - } - else{ - vxaverage_input=inputs->GetInput(VxAverageEnum); _assert_(vxaverage_input); - vyaverage_input=inputs->GetInput(VyAverageEnum); _assert_(vyaverage_input); - } - h=sqrt(2.*this->GetArea()); - - /*Start looping on the number of gaussian points:*/ - GaussTria* gauss=new GaussTria(2); - for(int ig=gauss->begin();igend();ig++){ - - gauss->GaussPoint(ig); - - GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss); - GetBMasstransport(B,&xyz_list[0][0],gauss); - GetBprimeMasstransport(Bprime,&xyz_list[0][0],gauss); - - vxaverage_input->GetInputValue(&vx,gauss); - vyaverage_input->GetInputValue(&vy,gauss); - vxaverage_input->GetInputDerivativeValue(&dvx[0],&xyz_list[0][0],gauss); - vyaverage_input->GetInputDerivativeValue(&dvy[0],&xyz_list[0][0],gauss); - - dvxdx=dvx[0]; - dvydy=dvy[1]; - D_scalar=gauss->weight*Jdet; - - D[0][0]=D_scalar*dvxdx; - D[0][1]=0.; - D[1][0]=0.; - D[1][1]=D_scalar*dvydy; - TripleMultiply(B,2,numnodes,1, - &D[0][0],2,2,0, - B,2,numnodes,0, - &Ke->values[0],1); - - D[0][0]=D_scalar*vx; - D[1][1]=D_scalar*vy; - TripleMultiply(B,2,numnodes,1, - &D[0][0],2,2,0, - Bprime,2,numnodes,0, - &Ke->values[0],1); - - if(stabilization==1){ - /*Streamline upwinding*/ - vel=sqrt(vx*vx+vy*vy); - D[0][0]=h/(2*vel)*vx*vx; - D[1][0]=h/(2*vel)*vy*vx; - D[0][1]=h/(2*vel)*vx*vy; - D[1][1]=h/(2*vel)*vy*vy; - } - else if(stabilization==2){ - /*SSA*/ - vxaverage_input->GetInputAverage(&vx); - vyaverage_input->GetInputAverage(&vy); - D[0][0]=h/2.0*fabs(vx); - D[0][1]=0.; - D[1][0]=0.; - D[1][1]=h/2.0*fabs(vy); - } - if(stabilization==1 || stabilization==2){ - D[0][0]=D_scalar*D[0][0]; - D[1][0]=D_scalar*D[1][0]; - D[0][1]=D_scalar*D[0][1]; - D[1][1]=D_scalar*D[1][1]; - TripleMultiply(Bprime,2,numnodes,1, - &D[0][0],2,2,0, - Bprime,2,numnodes,0, - &Ke->values[0],1); - } - } - - /*Clean up and return*/ - xDelete(basis); - xDelete(B); - xDelete(Bprime); - delete gauss; - return Ke; -} -/*}}}*/ -/*FUNCTION Tria::CreateKMatrixBalancethickness_DG {{{*/ -ElementMatrix* Tria::CreateKMatrixBalancethickness_DG(void){ - - /*Intermediaries*/ - IssmDouble vx,vy,D_scalar,Jdet; - IssmDouble xyz_list[NUMVERTICES][3]; - - /*Fetch number of nodes for this finite element*/ - int numnodes = this->NumberofNodes(); - - /*Initialize Element matrix and vectors*/ - ElementMatrix* Ke = new ElementMatrix(nodes,numnodes,this->parameters,NoneApproximationEnum); - IssmDouble* basis = xNew(numnodes); - IssmDouble* B = xNew(2*numnodes); - IssmDouble* Bprime = xNew(2*numnodes); - IssmDouble D[2][2]; - - /*Retrieve all inputs and parameters*/ - ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); - Input* vx_input=inputs->GetInput(VxEnum); _assert_(vx_input); - Input* vy_input=inputs->GetInput(VyEnum); _assert_(vy_input); - - /*Start looping on the number of gaussian points:*/ - GaussTria* gauss=new GaussTria(2); - for(int ig=gauss->begin();igend();ig++){ - - gauss->GaussPoint(ig); - - GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss); - /*WARNING: B and Bprime are inverted compared to usual masstransport!!!!*/ - GetBMasstransport(Bprime,&xyz_list[0][0],gauss); - GetBprimeMasstransport(B,&xyz_list[0][0],gauss); - - vx_input->GetInputValue(&vx,gauss); - vy_input->GetInputValue(&vy,gauss); - - D_scalar=-gauss->weight*Jdet; - D[0][0]=D_scalar*vx; - D[0][1]=0.; - D[1][0]=0.; - D[1][1]=D_scalar*vy; - - TripleMultiply(B,2,numnodes,1, - &D[0][0],2,2,0, - Bprime,2,numnodes,0, - &Ke->values[0],1); - } - - /*Clean up and return*/ - xDelete(basis); - xDelete(B); - xDelete(Bprime); - delete gauss; - return Ke; -} -/*}}}*/ -/*FUNCTION Tria::CreateKMatrixBalancevelocity{{{*/ -ElementMatrix* Tria::CreateKMatrixBalancevelocity(void){ - - /*Intermediaries */ - IssmDouble xyz_list[NUMVERTICES][3]; - IssmDouble dhdt_g,mb_g,ms_g,Jdet; - IssmDouble h,gamma,thickness; - IssmDouble hnx,hny,dhnx[2],dhny[2]; - - /*Fetch number of nodes and dof for this finite element*/ - int numnodes = this->NumberofNodes(); - - /*Initialize Element matrix and vectors*/ - ElementMatrix* Ke = new ElementMatrix(nodes,numnodes,this->parameters,NoneApproximationEnum); - IssmDouble* B = xNew(2*numnodes); - IssmDouble* basis = xNew(numnodes); - IssmDouble* dbasis = xNew(2*numnodes); - IssmDouble* HNx = xNew(numnodes); - IssmDouble* HNy = xNew(numnodes); - IssmDouble* H = xNew(numnodes); - IssmDouble* Nx = xNew(numnodes); - IssmDouble* Ny = xNew(numnodes); - - /*Retrieve all Inputs and parameters: */ - ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); - Input* H_input =inputs->GetInput(ThicknessEnum); _assert_(H_input); - h=sqrt(2.*this->GetArea()); - - /*Get vector N for all nodes and build HNx and HNy*/ - GetInputListOnNodes(Nx,SurfaceSlopeXEnum); - GetInputListOnNodes(Ny,SurfaceSlopeYEnum); - GetInputListOnNodes(H,ThicknessEnum); - for(int i=0;ibegin();igend();ig++){ - - gauss->GaussPoint(ig); - - H_input->GetInputValue(&thickness,gauss); - if(thickness<50.) thickness=50.; - TriaRef::GetInputDerivativeValue(&dhnx[0],HNx,&xyz_list[0][0],gauss); - TriaRef::GetInputDerivativeValue(&dhny[0],HNy,&xyz_list[0][0],gauss); - TriaRef::GetInputValue(&hnx,HNx,gauss); - TriaRef::GetInputValue(&hny,HNy,gauss); - - gamma=h/(2.*thickness+1.e-10); - - GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss); - GetNodalFunctions(basis,gauss); - GetNodalFunctionsDerivatives(dbasis,&xyz_list[0][0],gauss); - - for(int i=0;ivalues[i*numnodes+j] += gauss->weight*Jdet*( - (basis[i]+gamma*(basis[i]*(dhnx[0]+dhny[1]) + dbasis[0*numnodes+i]*hnx + dbasis[1*numnodes+i]*hny))* - (basis[j]*(dhnx[0]+dhny[1]) + dbasis[0*numnodes+j]*hnx + dbasis[1*numnodes+j]*hny) - ); - } - } - } - - /*Clean up and return*/ - xDelete(basis); - xDelete(dbasis); - xDelete(H); - xDelete(Nx); - xDelete(Ny); - xDelete(HNx); - xDelete(HNy); - xDelete(B); - delete gauss; - return Ke; -} -/*}}}*/ -/*FUNCTION Tria::CreateKMatrixSmoothedSlope {{{*/ -ElementMatrix* Tria::CreateKMatrixSmoothedSlope(void){ - - /* Intermediaries */ - IssmDouble D_scalar,Jdet,thickness; - IssmDouble xyz_list[NUMVERTICES][3]; - IssmDouble l=8.; - - /*Fetch number of nodes and dof for this finite element*/ - int numnodes = this->NumberofNodes(); - - /*Initialize Element matrix and vectors*/ - ElementMatrix* Ke = new ElementMatrix(nodes,numnodes,this->parameters); - IssmDouble* dbasis = xNew(2*numnodes); - IssmDouble* basis = xNew(numnodes); - - /*Retrieve all inputs and parameters*/ - ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); - Input* thickness_input=inputs->GetInput(ThicknessEnum); _assert_(thickness_input); - - /* Start looping on the number of gaussian points: */ - GaussTria* gauss=new GaussTria(2); - for(int ig=gauss->begin();igend();ig++){ - - gauss->GaussPoint(ig); - - GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss); - thickness_input->GetInputValue(&thickness,gauss); - if(thickness<50.) thickness=50.; - - GetNodalFunctions(basis,gauss); - GetNodalFunctionsDerivatives(dbasis,&xyz_list[0][0],gauss); - - for(int i=0;ivalues[i*numnodes+j] += gauss->weight*Jdet*( - basis[i]*basis[j] - +(l*thickness)*(l*thickness)*(dbasis[0*numnodes+i]*dbasis[0*numnodes+j] + dbasis[1*numnodes+i]*dbasis[1*numnodes+j]) - ); - } - } - } - - /*Clean up and return*/ - delete gauss; - xDelete(dbasis); - xDelete(basis); - return Ke; } /*}}}*/ Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tria.h (revision 16992) +++ /issm/trunk-jpl/src/c/classes/Elements/Tria.h (revision 16993) @@ -68,5 +68,4 @@ void SetwiseNodeConnectivity(int* d_nz,int* o_nz,Node* node,bool* flags,int* flagsindices,int set1_enum,int set2_enum); void CreateDVector(Vector* df); - void CreateJacobianMatrix(Matrix* Jff); void Delta18oParameterization(void); void ElementSizes(IssmDouble* hx,IssmDouble* hy,IssmDouble* hz){_error_("not implemented yet");}; @@ -202,24 +201,4 @@ void AddInput(int input_enum, IssmDouble* values, int interpolation_enum); void AddMaterialInput(int input_enum, IssmDouble* values, int interpolation_enum); - ElementMatrix* CreateKMatrix(void); - ElementMatrix* CreateKMatrixBalancethickness(void); - ElementMatrix* CreateKMatrixBalancethickness_DG(void); - ElementMatrix* CreateKMatrixBalancethickness_CG(void); - ElementMatrix* CreateKMatrixBalancevelocity(void); - ElementMatrix* CreateKMatrixSmoothedSlope(void); - ElementMatrix* CreateKMatrixMelting(void); - ElementMatrix* CreateKMatrixMasstransport(void); - ElementMatrix* CreateKMatrixMasstransport_CG(void); - ElementMatrix* CreateKMatrixMasstransport_DG(void); - ElementMatrix* CreateKMatrixExtrusion(void); - ElementMatrix* CreateKMatrixExtrusionVolume(void); - ElementMatrix* CreateKMatrixExtrusionSurface(void); - ElementMatrix* CreateKMatrixExtrusionBed(void); - ElementMatrix* CreateKMatrixFreeSurfaceTop(void); - ElementMatrix* CreateKMatrixFreeSurfaceTop1D(void); - ElementMatrix* CreateKMatrixFreeSurfaceBase(void); - ElementMatrix* CreateKMatrixFreeSurfaceBase1D(void); - ElementMatrix* CreateMassMatrix(void); - ElementMatrix* CreateBasalMassMatrix(void); IssmDouble EnthalpyDiffusionParameter(IssmDouble enthalpy,IssmDouble pressure){_error_("not implemented");}; IssmDouble EnthalpyDiffusionParameterVolume(int numvertices,IssmDouble* enthalpy,IssmDouble* pressure){_error_("not implemented");}; @@ -283,19 +262,6 @@ #ifdef _HAVE_STRESSBALANCE_ - ElementMatrix* CreateKMatrixStressbalanceSSA(void); - ElementMatrix* CreateKMatrixStressbalanceSSAViscous(void); - ElementMatrix* CreateKMatrixStressbalanceSSAFriction(void); - ElementMatrix* CreateKMatrixStressbalanceSIA(void); - ElementMatrix* CreateKMatrixStressbalanceFS(void); - ElementMatrix* CreateKMatrixStressbalanceFSViscous(void); - ElementMatrix* CreateKMatrixStressbalanceFSFriction(void); void PVectorGLSstabilization(ElementVector* pe); - ElementMatrix* CreateJacobianStressbalanceSSA(void); IssmDouble GetYcoord(GaussTria* gauss); - #endif - - #ifdef _HAVE_CONTROL_ - ElementMatrix* CreateKMatrixAdjointBalancethickness(void); - ElementMatrix* CreateKMatrixAdjointSSA(void); #endif @@ -309,7 +275,4 @@ #ifdef _HAVE_HYDROLOGY_ - ElementMatrix* CreateKMatrixHydrologyShreve(void); - ElementMatrix* CreateKMatrixHydrologyDCInefficient(void); - ElementMatrix* CreateKMatrixHydrologyDCEfficient(void); ElementMatrix* CreateEPLDomainMassMatrix(void); void CreateHydrologyWaterVelocityInput(void); @@ -324,6 +287,4 @@ #ifdef _HAVE_DAMAGE_ - ElementMatrix* CreateKMatrixDamageEvolution(void); - ElementMatrix* CreateKMatrixDamageEvolution_CG(void); void DamageEvolutionF(IssmDouble* flist); #endif Index: /issm/trunk-jpl/src/c/modules/CreateJacobianMatrixx/CreateJacobianMatrixx.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/CreateJacobianMatrixx/CreateJacobianMatrixx.cpp (revision 16992) +++ /issm/trunk-jpl/src/c/modules/CreateJacobianMatrixx/CreateJacobianMatrixx.cpp (revision 16993) @@ -10,6 +10,6 @@ void CreateJacobianMatrixx(Matrix** pJff,FemModel* femmodel,IssmDouble kmax){ - int i,connectivity; - int configuration_type; + int i; + int configuration_type,analysisenum; Element *element = NULL; Load *load = NULL; @@ -21,5 +21,6 @@ /*Recover some parameters*/ femmodel->parameters->FindParam(&configuration_type,ConfigurationTypeEnum); - femmodel->parameters->FindParam(&connectivity,MeshAverageVertexConnectivityEnum); + femmodel->parameters->FindParam(&analysisenum,AnalysisTypeEnum); + Analysis* analysis = EnumToAnalysis(analysisenum); /*Initialize Jacobian Matrix*/ @@ -29,5 +30,7 @@ for(i=0;ielements->Size();i++){ element=dynamic_cast(femmodel->elements->GetObjectByOffset(i)); - element->CreateJacobianMatrix(Jff); + ElementMatrix* Je = analysis->CreateJacobianMatrix(element); + if(Je) Je->AddToGlobal(Jff); + delete Je; } for (i=0;iloads->Size();i++){ @@ -39,4 +42,5 @@ /*Assign output pointer*/ + delete analysis; *pJff=Jff; Index: /issm/trunk-jpl/src/c/modules/CreateJacobianMatrixx/CreateJacobianMatrixx.h =================================================================== --- /issm/trunk-jpl/src/c/modules/CreateJacobianMatrixx/CreateJacobianMatrixx.h (revision 16992) +++ /issm/trunk-jpl/src/c/modules/CreateJacobianMatrixx/CreateJacobianMatrixx.h (revision 16993) @@ -6,4 +6,5 @@ #include "../../classes/classes.h" +#include "../../analyses/analyses.h" /* local prototypes: */ Index: /issm/trunk-jpl/src/c/modules/SystemMatricesx/SystemMatricesx.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/SystemMatricesx/SystemMatricesx.cpp (revision 16992) +++ /issm/trunk-jpl/src/c/modules/SystemMatricesx/SystemMatricesx.cpp (revision 16993) @@ -41,6 +41,5 @@ for (i=0;ielements->Size();i++){ element=dynamic_cast(femmodel->elements->GetObjectByOffset(i)); - ElementMatrix* Ke = element->CreateKMatrix(); - //ElementVector* pe = element->CreatePVector(); + ElementMatrix* Ke = analysis->CreateKMatrix(element); ElementVector* pe = analysis->CreatePVector(element); element->ReduceMatrices(Ke,pe); @@ -73,7 +72,5 @@ for (i=0;ielements->Size();i++){ element=dynamic_cast(femmodel->elements->GetObjectByOffset(i)); - ElementMatrix* Ke = element->CreateKMatrix(); - //ElementMatrix* Ke = analysis->CreateKMatrix(element); - //ElementVector* pe = element->CreatePVector(); + ElementMatrix* Ke = analysis->CreateKMatrix(element); ElementVector* pe = analysis->CreatePVector(element); element->ReduceMatrices(Ke,pe);