Index: /issm/trunk-jpl/src/c/classes/Elements/Element.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Element.h (revision 15059) +++ /issm/trunk-jpl/src/c/classes/Elements/Element.h (revision 15060) @@ -41,4 +41,5 @@ virtual void CreatePVector(Vector* pf)=0; virtual void CreateJacobianMatrix(Matrix* Jff)=0; + virtual void BasisIntegral(Vector* basisg)=0; virtual void GetSolutionFromInputs(Vector* solution)=0; virtual int GetNodeIndex(Node* node)=0; @@ -135,5 +136,4 @@ #ifdef _HAVE_HYDROLOGY_ virtual void GetHydrologyDCInefficientHmax(IssmDouble* ph_max, Node* innode)=0; - virtual void BasisIntegral(Vector* basisg)=0; virtual void GetHydrologyTransfer(Vector* transfer)=0; #endif Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp (revision 15059) +++ /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp (revision 15060) @@ -802,4 +802,40 @@ /*return output*/ return penta; +} +/*}}}*/ +/*FUNCTION Penta::BasisIntegral {{{*/ +void Penta::BasisIntegral(Vector* basisg){ + + /*Constants*/ + const int numdof=NDOF1*NUMVERTICES; + + /*Intermediaries */ + IssmDouble Jdet; + IssmDouble xyz_list[NUMVERTICES][3]; + IssmDouble basis[numdof]; + IssmDouble basisint[numdof]={0.}; + int *doflist=NULL; + GaussPenta* gauss=NULL; + + GetDofList(&doflist,NoneApproximationEnum,GsetEnum); + 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); + GetNodalFunctionsP1(&basis[0], gauss); + + for(int i=0;iweight*basis[i]; + } + + basisg->SetValues(numdof,doflist,&basisint[0],ADD_VAL); + + /*Clean up and return*/ + delete gauss; + xDelete(doflist); } /*}}}*/ @@ -2073,35 +2109,55 @@ void Penta::InputUpdateFromVector(IssmDouble* vector, int name, int type){ + const int numdof = NDOF1 *NUMVERTICES; + int *doflist = NULL; + IssmDouble values[numdof]; + /*Check that name is an element input*/ if (!IsInput(name)) return; - /*Penta update B in InputUpdateFromSolutionThermal, so don't look for B update here.*/ - switch(type){ - - case VertexEnum: - { - - /*New PentaVertexInpu*/ - IssmDouble values[6]; - - /*Get values on the 6 vertices*/ - for (int i=0;i<6;i++){ - values[i]=vector[this->vertices[i]->Pid()]; - } - - /*update input*/ - if (name==MaterialsRheologyZEnum || name==MaterialsRheologyZbarEnum){ - material->inputs->AddInput(new PentaP1Input(name,values)); - } - else{ - this->inputs->AddInput(new PentaP1Input(name,values)); - } - return; - break; + case VertexPIdEnum: + for (int i=0;ivertices[i]->Pid()]; } + /*update input*/ + if (name==MaterialsRheologyZEnum || name==MaterialsRheologyZbarEnum){ + material->inputs->AddInput(new PentaP1Input(name,values)); + } + else{ + this->inputs->AddInput(new PentaP1Input(name,values)); + } + return; + + case VertexSIdEnum: + for (int i=0;ivertices[i]->Sid()]; + } + /*update input*/ + if (name==MaterialsRheologyZEnum || name==MaterialsRheologyZbarEnum){ + material->inputs->AddInput(new PentaP1Input(name,values)); + } + else{ + this->inputs->AddInput(new PentaP1Input(name,values)); + } + return; + + case NodesEnum: + /*Get dof list: */ + GetDofList(&doflist,NoneApproximationEnum,GsetEnum); + + /*Use the dof list to index into the vector: */ + for(int i=0;i(values[i])) _error_("NaN found in solution vector"); + } + /*Add input to the element: */ + this->inputs->AddInput(new PentaP1Input(name,values)); + + /*Free ressources:*/ + xDelete(doflist); + return; default: - _error_("type " << type << " (" << EnumToStringx(type) << ") not implemented yet"); } @@ -2148,4 +2204,5 @@ name==InversionVyObsEnum || name==InversionVzObsEnum || + name==BasisIntegralEnum || name==TemperatureEnum || name==EnthalpyEnum || @@ -4514,12 +4571,12 @@ const int numdof=NDOF1*NUMVERTICES; - bool converged; - int i,rheology_law; - IssmDouble xyz_list[NUMVERTICES][3]; - IssmDouble values[numdof]; - IssmDouble B[numdof]; - IssmDouble B_average,s_average; - int* doflist=NULL; - //IssmDouble pressure[numdof]; + bool converged; + int i,rheology_law; + IssmDouble xyz_list[NUMVERTICES][3]; + IssmDouble values[numdof]; + IssmDouble B[numdof]; + IssmDouble B_average,s_average; + int *doflist = NULL; + IssmDouble pressure[numdof]; /*Get dof list: */ @@ -9318,11 +9375,4 @@ } /*}}}*/ -/*FUNCTION Penta::BasisIntegral {{{*/ -void Penta::BasisIntegral(Vector* basisg){ - _error_("Hydrological stuff not suported in Penta"); -} -/*}}}*/ - -/*}}}*/ /*FUNCTION Penta::GetHydrologyTransfer{{{*/ void Penta::GetHydrologyTransfer(Vector* transfer){ @@ -9330,4 +9380,3 @@ } /*}}}*/ - #endif Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Penta.h (revision 15059) +++ /issm/trunk-jpl/src/c/classes/Elements/Penta.h (revision 15060) @@ -265,4 +265,5 @@ void InputUpdateFromSolutionDiagnosticVert( IssmDouble* solutiong); void InputUpdateFromSolutionDiagnosticStokes( IssmDouble* solutiong); + void BasisIntegral(Vector* gbasis); void GetSolutionFromInputsDiagnosticHoriz(Vector* solutiong); void GetSolutionFromInputsDiagnosticHutter(Vector* solutiong); @@ -307,5 +308,4 @@ void CreateHydrologyWaterVelocityInput(void); void GetHydrologyDCInefficientHmax(IssmDouble* ph_max, Node* innode); - void BasisIntegral(Vector* gbasis); void GetHydrologyTransfer(Vector* transfer); #endif Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 15059) +++ /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 15060) @@ -1099,4 +1099,40 @@ _assert_(x2*y3 - y2*x3 + x1*y2 - y1*x2 + x3*y1 - y3*x1>0); return (x2*y3 - y2*x3 + x1*y2 - y1*x2 + x3*y1 - y3*x1)/2; +} +/*}}}*/ +/*FUNCTION Tria::BasisIntegral {{{*/ +void Tria::BasisIntegral(Vector* basisg){ + + /*Constants*/ + const int numdof=NDOF1*NUMVERTICES; + + /*Intermediaries */ + IssmDouble Jdet; + IssmDouble xyz_list[NUMVERTICES][3]; + IssmDouble basis[numdof]; + IssmDouble basisint[numdof]={0.}; + int *doflist=NULL; + GaussTria* gauss=NULL; + + GetDofList(&doflist,NoneApproximationEnum,GsetEnum); + 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); + GetNodalFunctions(&basis[0], gauss); + + for(int i=0;iweight*basis[i]; + } + + basisg->SetValues(numdof,doflist,&basisint[0],ADD_VAL); + + /*Clean up and return*/ + delete gauss; + xDelete(doflist); } /*}}}*/ @@ -1938,6 +1974,5 @@ switch(type){ case VertexPIdEnum: - /*Get values on the 3 vertices*/ - for (int i=0;i<3;i++){ + for (int i=0;ivertices[i]->Pid()]; } @@ -1952,6 +1987,5 @@ case VertexSIdEnum: - /*Get values on the 3 vertices*/ - for (int i=0;i<3;i++){ + for (int i=0;ivertices[i]->Sid()]; } @@ -6487,40 +6521,4 @@ } /*}}}*/ -/*FUNCTION Tria::BasisIntegral {{{*/ -void Tria::BasisIntegral(Vector* basisg){ - - /*Constants*/ - const int numdof=NDOF1*NUMVERTICES; - - /*Intermediaries */ - IssmDouble Jdet; - IssmDouble xyz_list[NUMVERTICES][3]; - IssmDouble basis[numdof]; - IssmDouble basisint[numdof]={0.}; - int *doflist=NULL; - GaussTria* gauss=NULL; - - GetDofList(&doflist,NoneApproximationEnum,GsetEnum); - 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); - GetNodalFunctions(&basis[0], gauss); - - for(int i=0;iweight*basis[i]; - } - - basisg->SetValues(numdof,doflist,&basisint[0],ADD_VAL); - - /*Clean up and return*/ - delete gauss; - xDelete(doflist); -} -/*}}}*/ /*FUNCTION Tria::GetHydrologyTransfer{{{*/ void Tria::GetHydrologyTransfer(Vector* transfer){ @@ -6582,5 +6580,4 @@ /*}}}*/ - #endif Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tria.h (revision 15059) +++ /issm/trunk-jpl/src/c/classes/Elements/Tria.h (revision 15060) @@ -198,4 +198,5 @@ ElementVector* CreatePVectorSlope(void); IssmDouble GetArea(void); + void BasisIntegral(Vector* gbasis); int GetElementType(void); void GetDofList(int** pdoflist,int approximation_enum,int setenum); @@ -255,5 +256,4 @@ void InputUpdateFromSolutionHydrologyDCEfficient(IssmDouble* solution); void GetHydrologyDCInefficientHmax(IssmDouble* ph_max, Node* innode); - void BasisIntegral(Vector* gbasis); void GetHydrologyTransfer(Vector* transfer); #endif