Index: /issm/trunk-jpl/src/c/analyses/L2ProjectionBaseAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/L2ProjectionBaseAnalysis.cpp (revision 16804) +++ /issm/trunk-jpl/src/c/analyses/L2ProjectionBaseAnalysis.cpp (revision 16805) @@ -59,5 +59,71 @@ }/*}}}*/ ElementVector* L2ProjectionBaseAnalysis::CreatePVector(Element* element){/*{{{*/ -_error_("not implemented yet"); + + /*Intermediaries*/ + int meshtype; + Element* basalelement; + + /*Get basal element*/ + element->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 vector*/ + ElementVector* pe = basalelement->NewElementVector(); + IssmDouble* basis = xNew(numnodes); + + /*Retrieve all inputs and parameters*/ + basalelement->GetVerticesCoordinates(&xyz_list); + basalelement->FindParam(&input_enum,InputToL2ProjectEnum); + switch(input_enum){ + case SurfaceSlopeXEnum: input2 = basalelement->GetInput(SurfaceEnum); _assert_(input2); break; + case SurfaceSlopeYEnum: input2 = basalelement->GetInput(SurfaceEnum); _assert_(input2); break; + case BedSlopeXEnum: input2 = basalelement->GetInput(BedEnum); _assert_(input2); break; + case BedSlopeYEnum: input2 = basalelement->GetInput(BedEnum); _assert_(input2); break; + default: input = element->GetInput(input_enum); + } + + /* 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->NodalFunctions(basis,gauss); + + if(input2) input2->GetInputDerivativeValue(&slopes[0],xyz_list,gauss); + switch(input_enum){ + case SurfaceSlopeXEnum: case BedSlopeXEnum: value = slopes[0]; break; + case SurfaceSlopeYEnum: case BedSlopeYEnum: value = slopes[1]; break; + default: input->GetInputValue(&value,gauss); + } + + for(int i=0;ivalues[i]+=Jdet*gauss->weight*value*basis[i]; + } + + /*Clean up and return*/ + xDelete(xyz_list); + xDelete(basis); + delete gauss; + if(meshtype!=Mesh2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;}; + return pe; }/*}}}*/ void L2ProjectionBaseAnalysis::GetSolutionFromInputs(Vector* solution,Element* element){/*{{{*/ Index: /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp (revision 16804) +++ /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp (revision 16805) @@ -814,7 +814,104 @@ case HOApproximationEnum: return CreatePVectorHO(element); + case FSApproximationEnum: + return CreatePVectorFS(element); + case NoneApproximationEnum: + return NULL; default: _error_("Approximation "<FindParam(&meshtype,MeshTypeEnum); + switch(meshtype){ + case Mesh2DverticalEnum: dim = 2; break; + case Mesh3DEnum: dim = 3; break; + default: _error_("mesh "<GetNumberOfNodesVelocity(); + int pnumnodes = element->GetNumberOfNodesPressure(); + + /*Prepare coordinate system list*/ + int* cs_list = xNew(vnumnodes+pnumnodes); + if(dim==2) for(i=0;iNewElementVector(FSvelocityEnum); + IssmDouble* vbasis = xNew(vnumnodes); + + /*Retrieve all inputs and parameters*/ + element->GetVerticesCoordinates(&xyz_list); + Input* loadingforcex_input=element->GetInput(LoadingforceXEnum); _assert_(loadingforcex_input); + Input* loadingforcey_input=element->GetInput(LoadingforceYEnum); _assert_(loadingforcey_input); + Input* loadingforcez_input=element->GetInput(LoadingforceZEnum); _assert_(loadingforcez_input); + IssmDouble rho_ice =element->GetMaterialParameter(MaterialsRhoIceEnum); + IssmDouble gravity =element->GetMaterialParameter(ConstantsGEnum); + + /* 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->NodalFunctionsVelocity(vbasis,gauss); + + loadingforcex_input->GetInputValue(&forcex,gauss); + loadingforcey_input->GetInputValue(&forcey,gauss); + if(dim==3) loadingforcez_input->GetInputValue(&forcez,gauss); + + for(i=0;ivalues[i*dim+0] += +rho_ice*forcex *Jdet*gauss->weight*vbasis[i]; + pe->values[i*dim+1] += +rho_ice*forcey *Jdet*gauss->weight*vbasis[i]; + if(dim==3){ + pe->values[i*dim+2] += +rho_ice*forcez*Jdet*gauss->weight*vbasis[i]; + pe->values[i*dim+2] += -rho_ice*gravity*Jdet*gauss->weight*vbasis[i]; + } + else{ + pe->values[i*dim+1] += -rho_ice*gravity*Jdet*gauss->weight*vbasis[i]; + } + } + } + + /*Transform coordinate system*/ + element->TransformLoadVectorCoord(pe,cs_list); + + /*Clean up and return*/ + delete gauss; + xDelete(cs_list); + xDelete(vbasis); + xDelete(xyz_list); + return pe; +}/*}}}*/ +ElementVector* StressbalanceAnalysis::CreatePVectorFSShelf(Element* element){/*{{{*/ + + return NULL; +}/*}}}*/ +ElementVector* StressbalanceAnalysis::CreatePVectorFSFront(Element* element){/*{{{*/ + + return NULL; }/*}}}*/ ElementVector* StressbalanceAnalysis::CreatePVectorHO(Element* element){/*{{{*/ @@ -916,5 +1013,5 @@ /*Initialize Element vector and vectors*/ - ElementVector* pe=element->NewElementVector(SSAApproximationEnum); + ElementVector* pe = element->NewElementVector(SSAApproximationEnum); IssmDouble* basis = xNew(numnodes); @@ -1146,10 +1243,6 @@ /*Prepare coordinate system list*/ int* cs_list = xNew(vnumnodes+pnumnodes); - if(dim==2){ - for(i=0;iEnum()==GaussPentaEnum); this->GetNodalFunctions(basis,(GaussPenta*)gauss); + +} +/*}}}*/ +/*FUNCTION Penta::NodalFunctionsVelocity{{{*/ +void Penta::NodalFunctionsVelocity(IssmDouble* basis, Gauss* gauss){ + + _assert_(gauss->Enum()==GaussPentaEnum); + this->GetNodalFunctionsVelocity(basis,(GaussPenta*)gauss); + +} +/*}}}*/ +/*FUNCTION Penta::NodalFunctionsPressure{{{*/ +void Penta::NodalFunctionsPressure(IssmDouble* basis, Gauss* gauss){ + + _assert_(gauss->Enum()==GaussPentaEnum); + this->GetNodalFunctionsPressure(basis,(GaussPenta*)gauss); } Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Penta.h (revision 16804) +++ /issm/trunk-jpl/src/c/classes/Elements/Penta.h (revision 16805) @@ -247,4 +247,6 @@ ElementVector* NewElementVector(int approximation_enum); void NodalFunctions(IssmDouble* basis,Gauss* gauss); + void NodalFunctionsVelocity(IssmDouble* basis,Gauss* gauss); + void NodalFunctionsPressure(IssmDouble* basis,Gauss* gauss); IssmDouble MinEdgeLength(IssmDouble xyz_list[6][3]); void SetClone(int* minranks); Index: /issm/trunk-jpl/src/c/classes/Elements/Seg.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Seg.h (revision 16804) +++ /issm/trunk-jpl/src/c/classes/Elements/Seg.h (revision 16805) @@ -111,4 +111,6 @@ void JacobianDeterminantBase(IssmDouble* pJdet,IssmDouble* xyz_list_base,Gauss* gauss){_error_("not implemented yet");}; void NodalFunctions(IssmDouble* basis,Gauss* gauss){_error_("not implemented yet");}; + void NodalFunctionsVelocity(IssmDouble* basis,Gauss* gauss){_error_("not implemented yet");}; + void NodalFunctionsPressure(IssmDouble* basis,Gauss* gauss){_error_("not implemented yet");}; bool NoIceInElement(){_error_("not implemented yet");}; int NumberofNodesVelocity(void){_error_("not implemented yet");}; Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 16804) +++ /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 16805) @@ -2207,4 +2207,20 @@ _assert_(gauss->Enum()==GaussTriaEnum); this->GetNodalFunctions(basis,(GaussTria*)gauss); + +} +/*}}}*/ +/*FUNCTION Tria::NodalFunctionsVelocity{{{*/ +void Tria::NodalFunctionsVelocity(IssmDouble* basis, Gauss* gauss){ + + _assert_(gauss->Enum()==GaussTriaEnum); + this->GetNodalFunctionsVelocity(basis,(GaussTria*)gauss); + +} +/*}}}*/ +/*FUNCTION Tria::NodalFunctionsPressure{{{*/ +void Tria::NodalFunctionsPressure(IssmDouble* basis, Gauss* gauss){ + + _assert_(gauss->Enum()==GaussTriaEnum); + this->GetNodalFunctionsPressure(basis,(GaussTria*)gauss); } Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tria.h (revision 16804) +++ /issm/trunk-jpl/src/c/classes/Elements/Tria.h (revision 16805) @@ -274,4 +274,6 @@ ElementVector* NewElementVector(int approximation_enum); void NodalFunctions(IssmDouble* basis,Gauss* gauss); + void NodalFunctionsVelocity(IssmDouble* basis,Gauss* gauss); + void NodalFunctionsPressure(IssmDouble* basis,Gauss* gauss); void SetClone(int* minranks); Seg* SpawnSeg(int index1,int index2); Index: /issm/trunk-jpl/src/c/modules/SystemMatricesx/SystemMatricesx.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/SystemMatricesx/SystemMatricesx.cpp (revision 16804) +++ /issm/trunk-jpl/src/c/modules/SystemMatricesx/SystemMatricesx.cpp (revision 16805) @@ -41,5 +41,10 @@ for (i=0;ielements->Size();i++){ element=dynamic_cast(femmodel->elements->GetObjectByOffset(i)); - element->CreateKMatrix(Kff_temp,NULL); + ElementMatrix* Ke = element->CreateKMatrix(); + ElementVector* pe = element->CreatePVector(); + element->ReduceMatrices(Ke,pe); + if(Ke) Ke->AddToGlobal(Kff_temp,NULL); + delete Ke; + delete pe; } @@ -68,7 +73,7 @@ element=dynamic_cast(femmodel->elements->GetObjectByOffset(i)); element=dynamic_cast(femmodel->elements->GetObjectByOffset(i)); - //ElementVector* pe = analysis->CreatePVector(element); + ElementVector* pe = analysis->CreatePVector(element); //ElementMatrix* Ke = analysis->CreateKMatrix(element); - ElementVector* pe = element->CreatePVector(); + //ElementVector* pe = element->CreatePVector(); ElementMatrix* Ke = element->CreateKMatrix(); element->ReduceMatrices(Ke,pe);