Index: /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp (revision 16826) +++ /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp (revision 16827) @@ -844,6 +844,6 @@ /*compute all stiffness matrices for this element*/ - ElementMatrix* Ke1=CreateKMatrixSSAViscous(element); - ElementMatrix* Ke2=CreateKMatrixSSAFriction(element); + ElementMatrix* Ke1=CreateKMatrixSSAViscous(basalelement); + ElementMatrix* Ke2=CreateKMatrixSSAFriction(basalelement); ElementMatrix* Ke =new ElementMatrix(Ke1,Ke2); Index: /issm/trunk-jpl/src/c/analyses/StressbalanceVerticalAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/StressbalanceVerticalAnalysis.cpp (revision 16826) +++ /issm/trunk-jpl/src/c/analyses/StressbalanceVerticalAnalysis.cpp (revision 16827) @@ -97,5 +97,101 @@ /*Finite Element Analysis*/ ElementMatrix* StressbalanceVerticalAnalysis::CreateKMatrix(Element* element){/*{{{*/ - _error_("not implemented yet"); + + /*compute all stiffness matrices for this element*/ + ElementMatrix* Ke1=CreateKMatrixVolume(element); + ElementMatrix* Ke2=CreateKMatrixSurface(element); + ElementMatrix* Ke =new ElementMatrix(Ke1,Ke2); + //Ke1->Echo(); + //Ke2->Echo(); + //Ke=element->CreateKMatrix(); + //_error_("S"); + + /*clean-up and return*/ + delete Ke1; + delete Ke2; + return Ke; + +}/*}}}*/ +ElementMatrix* StressbalanceVerticalAnalysis::CreateKMatrixVolume(Element* element){/*{{{*/ + + /*Intermediaries*/ + IssmDouble D,Jdet; + IssmDouble *xyz_list = NULL; + + /*Fetch number of nodes and dof for this finite element*/ + int numnodes = element->GetNumberOfNodes(); + + /*Initialize Element matrix and vectors*/ + ElementMatrix* Ke = element->NewElementMatrix(NoneApproximationEnum); + IssmDouble* B = xNew(numnodes); + IssmDouble* Bprime = xNew(numnodes); + + /*Retrieve all inputs and parameters*/ + element->GetVerticesCoordinates(&xyz_list); + + /* 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); + this->GetB(B,element,xyz_list,gauss); + this->GetBprime(Bprime,element,xyz_list,gauss); + D=gauss->weight*Jdet; + + TripleMultiply(B,1,numnodes,1, + &D,1,1,0, + Bprime,1,numnodes,0, + &Ke->values[0],1); + } + + /*Clean up and return*/ + delete gauss; + xDelete(xyz_list); + xDelete(Bprime); + xDelete(B); + return Ke; + +}/*}}}*/ +ElementMatrix* StressbalanceVerticalAnalysis::CreateKMatrixSurface(Element* element){/*{{{*/ + + + if(!element->IsOnSurface()) return NULL; + + /*Intermediaries*/ + IssmDouble D,Jdet,normal[3]; + IssmDouble *xyz_list = NULL; + + /*Fetch number of nodes and dof for this finite element*/ + int numnodes = element->GetNumberOfNodes(); + + /*Initialize Element matrix and vectors*/ + ElementMatrix* Ke = element->NewElementMatrix(NoneApproximationEnum); + IssmDouble* basis = xNew(numnodes); + + /*Retrieve all inputs and parameters*/ + element->GetVerticesCoordinatesTop(&xyz_list); + + /* Start looping on the number of gaussian points: */ + Gauss* gauss = element->NewGaussTop(2); + element->NormalTop(&normal[0],xyz_list); + for(int ig=gauss->begin();igend();ig++){ + gauss->GaussPoint(ig); + + element->JacobianDeterminantTop(&Jdet,xyz_list,gauss); + element->NodalFunctions(basis,gauss); + D = -gauss->weight*Jdet*normal[2]; + + TripleMultiply( basis,1,numnodes,1, + &D,1,1,0, + basis,1,numnodes,0, + &Ke->values[0],1); + } + + /*Clean up and return*/ + delete gauss; + xDelete(xyz_list); + xDelete(basis); + return Ke; }/*}}}*/ ElementVector* StressbalanceVerticalAnalysis::CreatePVector(Element* element){/*{{{*/ @@ -221,4 +317,28 @@ xDelete(xyz_list_base); return pe; +}/*}}}*/ +void StressbalanceVerticalAnalysis::GetB(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/ + /* Compute B matrix. B=[dh1/dz dh2/dz dh3/dz dh4/dz dh5/dz dh6/dz]; + where hi is the interpolation function for node i.*/ + + /*Fetch number of nodes for this finite element*/ + int numnodes = element->GetNumberOfNodes(); + + /*Get nodal functions derivatives*/ + IssmDouble* dbasis=xNew(3*numnodes); + element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss); + + /*Build B: */ + for(int i=0;i(dbasis); +}/*}}}*/ +void StressbalanceVerticalAnalysis::GetBprime(IssmDouble* Bprime,Element* element,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/ + + element->NodalFunctions(Bprime,gauss); + }/*}}}*/ void StressbalanceVerticalAnalysis::GetSolutionFromInputs(Vector* solution,Element* element){/*{{{*/ Index: /issm/trunk-jpl/src/c/analyses/StressbalanceVerticalAnalysis.h =================================================================== --- /issm/trunk-jpl/src/c/analyses/StressbalanceVerticalAnalysis.h (revision 16826) +++ /issm/trunk-jpl/src/c/analyses/StressbalanceVerticalAnalysis.h (revision 16827) @@ -22,7 +22,11 @@ /*Finite element Analysis*/ ElementMatrix* CreateKMatrix(Element* element); + ElementMatrix* CreateKMatrixVolume(Element* element); + ElementMatrix* CreateKMatrixSurface(Element* element); ElementVector* CreatePVector(Element* element); ElementVector* CreatePVectorVolume(Element* element); ElementVector* CreatePVectorBase(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 16826) +++ /issm/trunk-jpl/src/c/classes/Elements/Element.h (revision 16827) @@ -61,4 +61,5 @@ virtual void NodalFunctionsPressure(IssmDouble* basis,Gauss* gauss)=0; virtual void NormalBase(IssmDouble* normal,IssmDouble* xyz_list)=0; + virtual void NormalTop(IssmDouble* normal,IssmDouble* xyz_list)=0; virtual IssmDouble PureIceEnthalpy(IssmDouble pressure)=0; virtual void TransformLoadVectorCoord(ElementVector* pe,int transformenum)=0; @@ -80,4 +81,5 @@ virtual void JacobianDeterminant(IssmDouble* Jdet, IssmDouble* xyz_list,Gauss* gauss)=0; virtual void JacobianDeterminantBase(IssmDouble* Jdet,IssmDouble* xyz_list_base,Gauss* gauss)=0; + virtual void JacobianDeterminantTop(IssmDouble* Jdet,IssmDouble* xyz_list_base,Gauss* gauss)=0; virtual void GetSolutionFromInputsOneDof(Vector* solution,int solutionenum)=0; virtual int GetNodeIndex(Node* node)=0; @@ -93,4 +95,5 @@ virtual bool IsNodeOnShelfFromFlags(IssmDouble* flags)=0; virtual bool IsOnBed()=0; + virtual bool IsOnSurface()=0; virtual void GetInputListOnNodes(IssmDouble* pvalue,int enumtype)=0; virtual void GetInputListOnNodes(IssmDouble* pvalue,int enumtype,IssmDouble defaultvalue)=0; @@ -104,4 +107,5 @@ virtual void GetVerticesCoordinates(IssmDouble** xyz_list)=0; virtual void GetVerticesCoordinatesBase(IssmDouble** xyz_list)=0; + virtual void GetVerticesCoordinatesTop(IssmDouble** xyz_list)=0; virtual void GetMaterialInputValue(IssmDouble* pvalue,Node* node,int enumtype)=0; @@ -125,4 +129,5 @@ virtual Gauss* NewGauss(int order)=0; virtual Gauss* NewGaussBase(int order)=0; + virtual Gauss* NewGaussTop(int order)=0; virtual ElementVector* NewElementVector(int approximation_enum=NoneApproximationEnum)=0; virtual ElementMatrix* NewElementMatrix(int approximation_enum=NoneApproximationEnum)=0; Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp (revision 16826) +++ /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp (revision 16827) @@ -162,28 +162,4 @@ /*}}}*/ -/*FUNCTION Penta::NormalBase {{{*/ -void Penta::NormalBase(IssmDouble* bed_normal,IssmDouble* xyz_list){ - - int i; - IssmDouble v13[3],v23[3]; - IssmDouble normal[3]; - IssmDouble normal_norm; - - for (i=0;i<3;i++){ - v13[i]=xyz_list[0*3+i]-xyz_list[2*3+i]; - v23[i]=xyz_list[1*3+i]-xyz_list[2*3+i]; - } - - normal[0]=v13[1]*v23[2]-v13[2]*v23[1]; - normal[1]=v13[2]*v23[0]-v13[0]*v23[2]; - normal[2]=v13[0]*v23[1]-v13[1]*v23[0]; - normal_norm=sqrt( pow(normal[0],2)+pow(normal[1],2)+pow(normal[2],2) ); - - /*Bed normal is opposite to surface normal*/ - bed_normal[0]=-normal[0]/normal_norm; - bed_normal[1]=-normal[1]/normal_norm; - bed_normal[2]=-normal[2]/normal_norm; -} -/*}}}*/ /*FUNCTION Penta::BasalFrictionCreateInput {{{*/ void Penta::BasalFrictionCreateInput(void){ @@ -1466,4 +1442,14 @@ }/*}}}*/ +/*FUNCTION Penta::GetVerticesCoordinatesTop(IssmDouble** pxyz_list){{{*/ +void Penta::GetVerticesCoordinatesTop(IssmDouble** pxyz_list){ + + IssmDouble* xyz_list = xNew(NUMVERTICES2D*3); + ::GetVerticesCoordinates(xyz_list,&this->vertices[3],NUMVERTICES2D); + + /*Assign output pointer*/ + *pxyz_list = xyz_list; + +}/*}}}*/ /*FUNCTION Penta::GetMaterialInputValue(IssmDouble* pvalue,Node* node,int enumtype) {{{*/ void Penta::GetMaterialInputValue(IssmDouble* pvalue,Node* node,int enumtype){ @@ -2334,4 +2320,12 @@ inputs->GetInputValue(&onbed,MeshElementonbedEnum); return onbed; +} +/*}}}*/ +/*FUNCTION Penta::IsOnSurface{{{*/ +bool Penta::IsOnSurface(void){ + + bool onsurface; + inputs->GetInputValue(&onsurface,MeshElementonsurfaceEnum); + return onsurface; } /*}}}*/ @@ -2436,12 +2430,4 @@ } /*}}}*/ -/*FUNCTION Penta::IsOnSurface{{{*/ -bool Penta::IsOnSurface(void){ - - bool onsurface; - inputs->GetInputValue(&onsurface,MeshElementonsurfaceEnum); - return onsurface; -} -/*}}}*/ /*FUNCTION Penta::JacobianDeterminant{{{*/ void Penta::JacobianDeterminant(IssmDouble* pJdet,IssmDouble* xyz_list,Gauss* gauss){ @@ -2457,4 +2443,12 @@ _assert_(gauss->Enum()==GaussPentaEnum); this->GetTriaJacobianDeterminant(pJdet,xyz_list_base,(GaussPenta*)gauss); + +} +/*}}}*/ +/*FUNCTION Penta::JacobianDeterminantTop{{{*/ +void Penta::JacobianDeterminantTop(IssmDouble* pJdet,IssmDouble* xyz_list_top,Gauss* gauss){ + + _assert_(gauss->Enum()==GaussPentaEnum); + this->GetTriaJacobianDeterminant(pJdet,xyz_list_top,(GaussPenta*)gauss); } @@ -2543,4 +2537,9 @@ } /*}}}*/ +/*FUNCTION Penta::NewGaussTop(int order){{{*/ +Gauss* Penta::NewGaussTop(int order){ + return new GaussPenta(3,4,5,order); +} +/*}}}*/ /*FUNCTION Penta::NewElementVector{{{*/ ElementVector* Penta::NewElementVector(int approximation_enum){ @@ -2583,4 +2582,51 @@ this->GetNodalFunctionsPressure(basis,(GaussPenta*)gauss); +} +/*}}}*/ +/*FUNCTION Penta::NormalBase {{{*/ +void Penta::NormalBase(IssmDouble* bed_normal,IssmDouble* xyz_list){ + + int i; + IssmDouble v13[3],v23[3]; + IssmDouble normal[3]; + IssmDouble normal_norm; + + for (i=0;i<3;i++){ + v13[i]=xyz_list[0*3+i]-xyz_list[2*3+i]; + v23[i]=xyz_list[1*3+i]-xyz_list[2*3+i]; + } + + normal[0]=v13[1]*v23[2]-v13[2]*v23[1]; + normal[1]=v13[2]*v23[0]-v13[0]*v23[2]; + normal[2]=v13[0]*v23[1]-v13[1]*v23[0]; + normal_norm=sqrt( pow(normal[0],2)+pow(normal[1],2)+pow(normal[2],2) ); + + /*Bed normal is opposite to surface normal*/ + bed_normal[0]=-normal[0]/normal_norm; + bed_normal[1]=-normal[1]/normal_norm; + bed_normal[2]=-normal[2]/normal_norm; +} +/*}}}*/ +/*FUNCTION Penta::NormalTop {{{*/ +void Penta::NormalTop(IssmDouble* top_normal,IssmDouble* xyz_list){ + + int i; + IssmDouble v13[3],v23[3]; + IssmDouble normal[3]; + IssmDouble normal_norm; + + for (i=0;i<3;i++){ + v13[i]=xyz_list[0*3+i]-xyz_list[2*3+i]; + v23[i]=xyz_list[1*3+i]-xyz_list[2*3+i]; + } + + normal[0]=v13[1]*v23[2]-v13[2]*v23[1]; + normal[1]=v13[2]*v23[0]-v13[0]*v23[2]; + normal[2]=v13[0]*v23[1]-v13[1]*v23[0]; + normal_norm=sqrt(normal[0]*normal[0] + normal[1]*normal[1] + normal[2]*normal[2]); + + top_normal[0]=normal[0]/normal_norm; + top_normal[1]=normal[1]/normal_norm; + top_normal[2]=normal[2]/normal_norm; } /*}}}*/ @@ -2954,5 +3000,11 @@ _assert_(this->IsOnBed()); + + this->InputDepthAverageAtBase(MaterialsRheologyBEnum,MaterialsRheologyBbarEnum,MaterialsEnum); + this->InputDepthAverageAtBase(DamageDEnum,DamageDbarEnum,MaterialsEnum); Tria* tria=(Tria*)SpawnTria(0); //lower face is 0, upper face is 1. + this->material->inputs->DeleteInput(MaterialsRheologyBbarEnum); + this->material->inputs->DeleteInput(DamageDbarEnum); + return tria; } Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Penta.h (revision 16826) +++ /issm/trunk-jpl/src/c/classes/Elements/Penta.h (revision 16827) @@ -100,4 +100,5 @@ void GetVerticesCoordinates(IssmDouble** pxyz_list); void GetVerticesCoordinatesBase(IssmDouble** pxyz_list); + void GetVerticesCoordinatesTop(IssmDouble** pxyz_list); int Sid(); @@ -195,4 +196,5 @@ void AddMaterialInput(int input_enum, IssmDouble* values, int interpolation_enum); void NormalBase(IssmDouble* bed_normal, IssmDouble* xyz_list); + void NormalTop(IssmDouble* bed_normal, IssmDouble* xyz_list); ElementMatrix* CreateBasalMassMatrix(void); ElementMatrix* CreateKMatrix(void); @@ -244,8 +246,10 @@ void JacobianDeterminant(IssmDouble* Jdet, IssmDouble* xyz_list,Gauss* gauss); void JacobianDeterminantBase(IssmDouble* pJdet,IssmDouble* xyz_list_base,Gauss* gauss); + void JacobianDeterminantTop(IssmDouble* pJdet,IssmDouble* xyz_list_base,Gauss* gauss); bool NoIceInElement(void); Gauss* NewGauss(void); Gauss* NewGauss(int order); Gauss* NewGaussBase(int order); + Gauss* NewGaussTop(int order); ElementVector* NewElementVector(int approximation_enum); ElementMatrix* NewElementMatrix(int approximation_enum); Index: /issm/trunk-jpl/src/c/classes/Elements/Seg.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Seg.h (revision 16826) +++ /issm/trunk-jpl/src/c/classes/Elements/Seg.h (revision 16827) @@ -105,11 +105,14 @@ void GetVerticesCoordinates(IssmDouble** pxyz_list){_error_("not implemented yet");}; void GetVerticesCoordinatesBase(IssmDouble** pxyz_list){_error_("not implemented yet");}; + void GetVerticesCoordinatesTop(IssmDouble** pxyz_list){_error_("not implemented yet");}; int Sid(){_error_("not implemented yet");}; void InputChangeName(int input_enum, int enum_type_old){_error_("not implemented yet");}; bool IsOnBed(){_error_("not implemented yet");}; + bool IsOnSurface(){_error_("not implemented yet");}; bool IsFloating(){_error_("not implemented yet");}; bool IsNodeOnShelfFromFlags(IssmDouble* flags){_error_("not implemented yet");}; void JacobianDeterminant(IssmDouble* Jdet, IssmDouble* xyz_list,Gauss* gauss){_error_("not implemented yet");}; void JacobianDeterminantBase(IssmDouble* pJdet,IssmDouble* xyz_list_base,Gauss* gauss){_error_("not implemented yet");}; + void JacobianDeterminantTop(IssmDouble* pJdet,IssmDouble* xyz_list_base,Gauss* gauss){_error_("not implemented yet");}; IssmDouble MinEdgeLength(IssmDouble* xyz_list){_error_("not implemented yet");}; void NodalFunctions(IssmDouble* basis,Gauss* gauss){_error_("not implemented yet");}; @@ -119,4 +122,5 @@ bool NoIceInElement(){_error_("not implemented yet");}; void NormalBase(IssmDouble* normal,IssmDouble* xyz_list){_error_("not implemented yet");}; + void NormalTop(IssmDouble* normal,IssmDouble* xyz_list){_error_("not implemented yet");}; int NumberofNodesVelocity(void){_error_("not implemented yet");}; int NumberofNodesPressure(void){_error_("not implemented yet");}; @@ -137,4 +141,5 @@ Gauss* NewGauss(int order){_error_("not implemented yet");}; Gauss* NewGaussBase(int order){_error_("not implemented yet");}; + Gauss* NewGaussTop(int order){_error_("not implemented yet");}; ElementVector* NewElementVector(int approximation_enum){_error_("not implemented yet");}; ElementMatrix* NewElementMatrix(int approximation_enum){_error_("not implemented yet");}; Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 16826) +++ /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 16827) @@ -1980,4 +1980,18 @@ case Mesh2DverticalEnum: return HasEdgeOnBed(); + case Mesh2DhorizontalEnum: + return true; + default: _error_("mesh "<parameters->FindParam(&meshtype,MeshTypeEnum); + switch(meshtype){ + case Mesh2DverticalEnum: + return HasEdgeOnSurface(); case Mesh2DhorizontalEnum: return true; Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tria.h (revision 16826) +++ /issm/trunk-jpl/src/c/classes/Elements/Tria.h (revision 16827) @@ -95,4 +95,5 @@ int Sid(); bool IsOnBed(); + bool IsOnSurface(); bool HasEdgeOnBed(); bool HasNodeOnBed(); @@ -112,4 +113,5 @@ void GetVerticesCoordinates(IssmDouble** pxyz_list); void GetVerticesCoordinatesBase(IssmDouble** pxyz_list){_error_("not implemented yet");}; + void GetVerticesCoordinatesTop(IssmDouble** pxyz_list){_error_("not implemented yet");}; void InputCreate(IssmDouble* vector,IoModel* iomodel,int M,int N,int vector_type,int vector_enum,int code); void InputDepthAverageAtBase(int enum_type,int average_enum_type,int object_enum=MeshElementsEnum); @@ -253,4 +255,5 @@ IssmDouble GetGroundedPortion(IssmDouble* xyz_list); void NormalBase(IssmDouble* normal,IssmDouble* xyz_list); + void NormalTop(IssmDouble* normal,IssmDouble* xyz_list){_error_("not implemented yet");}; IssmDouble GetMaterialParameter(int enum_in); void GetZeroLevelsetCoordinates(IssmDouble* xyz_zero,IssmDouble xyz_list[3][3],int levelsetenum); @@ -273,8 +276,10 @@ void JacobianDeterminant(IssmDouble* pJdet, IssmDouble* xyz_list,Gauss* gauss); void JacobianDeterminantBase(IssmDouble* pJdet,IssmDouble* xyz_list_base,Gauss* gauss){_error_("not implemented yet");}; + void JacobianDeterminantTop(IssmDouble* pJdet,IssmDouble* xyz_list_base,Gauss* gauss){_error_("not implemented yet");}; IssmDouble MinEdgeLength(IssmDouble* xyz_list){_error_("not implemented yet");}; Gauss* NewGauss(void); Gauss* NewGauss(int order); Gauss* NewGaussBase(int order){_error_("not implemented yet");}; + Gauss* NewGaussTop(int order){_error_("not implemented yet");}; ElementVector* NewElementVector(int approximation_enum); ElementMatrix* NewElementMatrix(int approximation_enum);