Index: /issm/trunk-jpl/src/c/analyses/BalancethicknessAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/BalancethicknessAnalysis.cpp (revision 16831) +++ /issm/trunk-jpl/src/c/analyses/BalancethicknessAnalysis.cpp (revision 16832) @@ -136,5 +136,5 @@ /*Intermediaries */ IssmDouble dhdt,mb,ms,Jdet; - IssmDouble* xyz_list; + IssmDouble* xyz_list = NULL; /*Fetch number of nodes and dof for this finite element*/ @@ -177,5 +177,5 @@ /*Intermediaries */ IssmDouble dhdt,mb,ms,Jdet; - IssmDouble* xyz_list; + IssmDouble* xyz_list = NULL; /*Fetch number of nodes and dof for this finite element*/ Index: /issm/trunk-jpl/src/c/analyses/MasstransportAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/MasstransportAnalysis.cpp (revision 16831) +++ /issm/trunk-jpl/src/c/analyses/MasstransportAnalysis.cpp (revision 16832) @@ -213,5 +213,208 @@ /*Finite Element Analysis*/ ElementMatrix* MasstransportAnalysis::CreateKMatrix(Element* element){/*{{{*/ - _error_("not implemented yet"); + + if(!element->IsOnBed()) return NULL; + Element* basalelement = element->SpawnBasalElement(); + + int meshtype; + element->FindParam(&meshtype,MeshTypeEnum); + + switch(element->FiniteElement()){ + case P1Enum: case P2Enum: + return CreateKMatrixCG(basalelement); + case P1DGEnum: + return CreateKMatrixDG(basalelement); + default: + _error_("Element type " << EnumToStringx(element->FiniteElement()) << " not supported yet"); + } +}/*}}}*/ +ElementMatrix* MasstransportAnalysis::CreateKMatrixCG(Element* element){/*{{{*/ + + /*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 = NULL; + + /*Fetch number of nodes and dof for this finite element*/ + int numnodes = element->GetNumberOfNodes(); + + /*Initialize Element vector and other vectors*/ + ElementMatrix* Ke = element->NewElementMatrix(); + IssmDouble* basis = xNew(numnodes); + IssmDouble* B = xNew(2*numnodes); + IssmDouble* Bprime = xNew(2*numnodes); + IssmDouble D[2][2]; + + /*Retrieve all inputs and parameters*/ + element->GetVerticesCoordinates(&xyz_list); + element->FindParam(&dt,TimesteppingTimeStepEnum); + element->FindParam(&meshtype,MeshTypeEnum); + element->FindParam(&stabilization,MasstransportStabilizationEnum); + Input* vxaverage_input=NULL; + Input* vyaverage_input=NULL; + if(meshtype==Mesh2DhorizontalEnum){ + vxaverage_input=element->GetInput(VxEnum); _assert_(vxaverage_input); + vyaverage_input=element->GetInput(VyEnum); _assert_(vyaverage_input); + } + else{ + vxaverage_input=element->GetInput(VxAverageEnum); _assert_(vxaverage_input); + vyaverage_input=element->GetInput(VyAverageEnum); _assert_(vyaverage_input); + } + h = element->CharacteristicLength(); + + /* Start looping on the number of gaussian points: */ + Gauss* gauss=element->NewGauss(2); + for(int ig=gauss->begin();igend();ig++){ + gauss->GaussPoint(ig); + + element->JacobianDeterminant(&Jdet,xyz_list,gauss); + element->NodalFunctions(basis,gauss); + + vxaverage_input->GetInputValue(&vx,gauss); + vyaverage_input->GetInputValue(&vy,gauss); + vxaverage_input->GetInputDerivativeValue(&dvx[0],xyz_list,gauss); + vyaverage_input->GetInputDerivativeValue(&dvy[0],xyz_list,gauss); + + D_scalar=gauss->weight*Jdet; + TripleMultiply(basis,1,numnodes,1, + &D_scalar,1,1,0, + basis,1,numnodes,0, + &Ke->values[0],1); + + GetB(B,element,xyz_list,gauss); + GetBprime(Bprime,element,xyz_list,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(xyz_list); + xDelete(basis); + xDelete(B); + xDelete(Bprime); + delete gauss; + return Ke; +}/*}}}*/ +ElementMatrix* MasstransportAnalysis::CreateKMatrixDG(Element* element){/*{{{*/ + + /*Intermediaries */ + int meshtype; + IssmDouble Jdet,D_scalar,dt,vx,vy; + IssmDouble* xyz_list = NULL; + + /*Fetch number of nodes and dof for this finite element*/ + int numnodes = element->GetNumberOfNodes(); + + /*Initialize Element vector and other vectors*/ + ElementMatrix* Ke = element->NewElementMatrix(); + IssmDouble* basis = xNew(numnodes); + IssmDouble* B = xNew(2*numnodes); + IssmDouble* Bprime = xNew(2*numnodes); + IssmDouble D[2][2]; + + /*Retrieve all inputs and parameters*/ + element->GetVerticesCoordinates(&xyz_list); + element->FindParam(&dt,TimesteppingTimeStepEnum); + element->FindParam(&meshtype,MeshTypeEnum); + Input* vxaverage_input=NULL; + Input* vyaverage_input=NULL; + if(meshtype==Mesh2DhorizontalEnum){ + vxaverage_input=element->GetInput(VxEnum); _assert_(vxaverage_input); + vyaverage_input=element->GetInput(VyEnum); _assert_(vyaverage_input); + } + else{ + vxaverage_input=element->GetInput(VxAverageEnum); _assert_(vxaverage_input); + vyaverage_input=element->GetInput(VyAverageEnum); _assert_(vyaverage_input); + } + + /* Start looping on the number of gaussian points: */ + Gauss* gauss=element->NewGauss(2); + for(int ig=gauss->begin();igend();ig++){ + gauss->GaussPoint(ig); + + element->JacobianDeterminant(&Jdet,xyz_list,gauss); + element->NodalFunctions(basis,gauss); + + vxaverage_input->GetInputValue(&vx,gauss); + vyaverage_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); + + /*WARNING: B and Bprime are inverted compared to CG*/ + GetB(Bprime,element,xyz_list,gauss); + GetBprime(B,element,xyz_list,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(xyz_list); + xDelete(basis); + xDelete(B); + xDelete(Bprime); + delete gauss; + return Ke; }/*}}}*/ ElementVector* MasstransportAnalysis::CreatePVector(Element* element){/*{{{*/ @@ -235,5 +438,5 @@ IssmDouble Jdet,dt; IssmDouble ms,mb,mb_correction=0.,thickness; - IssmDouble* xyz_list; + IssmDouble* xyz_list = NULL; /*Fetch number of nodes and dof for this finite element*/ @@ -281,5 +484,5 @@ IssmDouble Jdet,dt; IssmDouble ms,mb,mb_correction=0.,thickness; - IssmDouble* xyz_list; + IssmDouble* xyz_list = NULL; /*Fetch number of nodes and dof for this finite element*/ @@ -321,4 +524,59 @@ delete gauss; return pe; +}/*}}}*/ +void MasstransportAnalysis::GetB(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/ + /*Compute B matrix. B=[B1 B2 B3] where Bi is of size 3*NDOF2. + * For node i, Bi can be expressed in the actual coordinate system + * by: + * Bi=[ N ] + * [ N ] + * where N is the finiteelement function for node i. + * + * We assume B_prog has been allocated already, of size: 2x(NDOF1*numnodes) + */ + + /*Fetch number of nodes for this finite element*/ + int numnodes = element->GetNumberOfNodes(); + + /*Get nodal functions*/ + IssmDouble* basis=xNew(numnodes); + element->NodalFunctions(basis,gauss); + + /*Build B: */ + for(int i=0;i(basis); +}/*}}}*/ +void MasstransportAnalysis::GetBprime(IssmDouble* Bprime,Element* element,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/ + /*Compute B' matrix. B'=[B1' B2' B3'] where Bi' is of size 3*NDOF2. + * For node i, Bi' can be expressed in the actual coordinate system + * by: + * Bi_prime=[ dN/dx ] + * [ dN/dy ] + * where N is the finiteelement function for node i. + * + * We assume B' has been allocated already, of size: 3x(NDOF2*numnodes) + */ + + /*Fetch number of nodes for this finite element*/ + int numnodes = element->GetNumberOfNodes(); + + /*Get nodal functions derivatives*/ + IssmDouble* dbasis=xNew(2*numnodes); + element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss); + + /*Build B': */ + for(int i=0;i(dbasis); + }/*}}}*/ void MasstransportAnalysis::GetSolutionFromInputs(Vector* solution,Element* element){/*{{{*/ Index: /issm/trunk-jpl/src/c/analyses/MasstransportAnalysis.h =================================================================== --- /issm/trunk-jpl/src/c/analyses/MasstransportAnalysis.h (revision 16831) +++ /issm/trunk-jpl/src/c/analyses/MasstransportAnalysis.h (revision 16832) @@ -22,7 +22,11 @@ /*Finite element Analysis*/ ElementMatrix* CreateKMatrix(Element* element); + ElementMatrix* CreateKMatrixCG(Element* element); + ElementMatrix* CreateKMatrixDG(Element* element); ElementVector* CreatePVector(Element* element); ElementVector* CreatePVectorCG(Element* element); ElementVector* CreatePVectorDG(Element* element); + void GetB(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss); + void GetBprime(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss); void GetSolutionFromInputs(Vector* solution,Element* element); void InputUpdateFromSolution(IssmDouble* solution,Element* element); Index: /issm/trunk-jpl/src/c/classes/Elements/Element.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Element.h (revision 16831) +++ /issm/trunk-jpl/src/c/classes/Elements/Element.h (revision 16832) @@ -37,4 +37,5 @@ virtual ~Element(){}; + virtual IssmDouble CharacteristicLength(void)=0; virtual void Configure(Elements* elements,Loads* loads,Nodes* nodes,Vertices* vertices,Materials* materials,Parameters* parameters)=0; virtual void SetCurrentConfiguration(Elements* elements,Loads* loads,Nodes* nodes,Materials* materials,Parameters* parameters)=0; Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp (revision 16831) +++ /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp (revision 16832) @@ -3011,7 +3011,11 @@ this->InputDepthAverageAtBase(MaterialsRheologyBEnum,MaterialsRheologyBbarEnum,MaterialsEnum); this->InputDepthAverageAtBase(DamageDEnum,DamageDbarEnum,MaterialsEnum); + if(this->inputs->GetInput(VxEnum)) this->InputDepthAverageAtBase(VxEnum,VxAverageEnum); + if(this->inputs->GetInput(VyEnum)) this->InputDepthAverageAtBase(VyEnum,VyAverageEnum); Tria* tria=(Tria*)SpawnTria(0); //lower face is 0, upper face is 1. this->material->inputs->DeleteInput(MaterialsRheologyBbarEnum); this->material->inputs->DeleteInput(DamageDbarEnum); + this->inputs->DeleteInput(VxAverageEnum); + this->inputs->DeleteInput(VyAverageEnum); return tria; Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Penta.h (revision 16831) +++ /issm/trunk-jpl/src/c/classes/Elements/Penta.h (revision 16832) @@ -69,4 +69,5 @@ /*Element virtual functions definitions: {{{*/ void BasalFrictionCreateInput(void); + IssmDouble CharacteristicLength(void){_error_("not implemented yet");}; void ComputeBasalStress(Vector* sigma_b); void ComputeStrainRate(Vector* eps); Index: /issm/trunk-jpl/src/c/classes/Elements/Seg.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Seg.h (revision 16831) +++ /issm/trunk-jpl/src/c/classes/Elements/Seg.h (revision 16832) @@ -70,4 +70,5 @@ void AddInput(int input_enum, IssmDouble* values, int interpolation_enum){_error_("not implemented yet");}; void AddMaterialInput(int input_enum, IssmDouble* values, int interpolation_enum){_error_("not implemented yet");}; + IssmDouble CharacteristicLength(void){_error_("not implemented yet");}; void ComputeBasalStress(Vector* sigma_b){_error_("not implemented yet");}; void ComputeStrainRate(Vector* eps){_error_("not implemented yet");}; Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 16831) +++ /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 16832) @@ -196,4 +196,10 @@ _assert_(this->material); this->material->inputs->AddInput(new TriaInput(input_enum,values,interpolation_enum)); +} +/*}}}*/ +/*FUNCTION Tria::AddInput{{{*/ +IssmDouble Tria::CharacteristicLength(void){ + + return sqrt(2*this->GetArea()); } /*}}}*/ Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tria.h (revision 16831) +++ /issm/trunk-jpl/src/c/classes/Elements/Tria.h (revision 16832) @@ -67,4 +67,5 @@ /*}}}*/ /*Element virtual functions definitions: {{{*/ + IssmDouble CharacteristicLength(void); void ComputeBasalStress(Vector* sigma_b); void ComputeStrainRate(Vector* eps);