Index: /issm/trunk/src/c/objects/Elements/Element.h =================================================================== --- /issm/trunk/src/c/objects/Elements/Element.h (revision 5634) +++ /issm/trunk/src/c/objects/Elements/Element.h (revision 5635) @@ -32,5 +32,4 @@ virtual void GetSolutionFromInputs(Vec solution)=0; virtual void GetNodes(void** nodes)=0; - virtual void* GetMatPar()=0; virtual bool GetShelf()=0; virtual bool GetOnBed()=0; Index: /issm/trunk/src/c/objects/Elements/Penta.cpp =================================================================== --- /issm/trunk/src/c/objects/Elements/Penta.cpp (revision 5634) +++ /issm/trunk/src/c/objects/Elements/Penta.cpp (revision 5635) @@ -941,10 +941,4 @@ this->results=new Results(); -} -/*}}}*/ -/*FUNCTION Penta::GetMatPar {{{1*/ -void* Penta::GetMatPar(){ - - return matpar; } /*}}}*/ @@ -4860,5 +4854,5 @@ /*clean-up*/ - delete gauss; + //delete gauss; } /*}}}*/ Index: /issm/trunk/src/c/objects/Elements/Penta.h =================================================================== --- /issm/trunk/src/c/objects/Elements/Penta.h (revision 5634) +++ /issm/trunk/src/c/objects/Elements/Penta.h (revision 5635) @@ -79,5 +79,4 @@ void CreatePVector(Vec pg); void DeleteResults(void); - void* GetMatPar(); void GetNodes(void** nodes); bool GetOnBed(); Index: /issm/trunk/src/c/objects/Elements/Tria.cpp =================================================================== --- /issm/trunk/src/c/objects/Elements/Tria.cpp (revision 5634) +++ /issm/trunk/src/c/objects/Elements/Tria.cpp (revision 5635) @@ -572,14 +572,5 @@ /*FUNCTION Tria::ComputeBasalStress {{{1*/ void Tria::ComputeBasalStress(Vec eps){ - - int i; - const int numvertices=3; - int doflist[numvertices]; - double value; - - /*plug local pressure values into global pressure vector: */ ISSMERROR("Not Implemented yet"); - //VecSetValues(eps,1,2,(const double*)&value,INSERT_VALUES); - } /*}}}*/ @@ -587,10 +578,9 @@ void Tria::ComputePressure(Vec pg){ - int i; - const int numvertices=3; - int doflist[numvertices]; - double pressure[numvertices]; - double thickness[numvertices]; - double rho_ice,g; + const int numvertices= 3; + int doflist[numvertices]; + double pressure[numvertices]; + double thickness[numvertices]; + double rho_ice,g; /*Get dof list on which we will plug the pressure values: */ @@ -602,5 +592,5 @@ GetParameterListOnVertices(&thickness[0],ThicknessEnum); - for(i=0;iGetParameterValue(&onwater,ElementOnWaterEnum); + + /*Intermediaries*/ + int i,ig; + bool onwater; + double thickness,thicknessobs,weight; + double Jdet; + double Jelem = 0; + double xyz_list[numvertices][3]; + GaussTria *gauss = NULL; /*If on water, return 0: */ @@ -1825,37 +1776,29 @@ /*Retrieve all inputs we will be needing: */ - thickness_input =inputs->GetInput(ThicknessEnum); ISSMASSERT(thickness_input); - thicknessobs_input=inputs->GetInput(ThicknessObsEnum);ISSMASSERT(thicknessobs_input); - weights_input =inputs->GetInput(WeightsEnum);ISSMASSERT(weights_input); - - /* Get gaussian points and weights (make this a statically initialized list of points? fstd): */ - GaussLegendreTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2); + inputs->GetParameterValue(&onwater,ElementOnWaterEnum); + Input* thickness_input =inputs->GetInput(ThicknessEnum); ISSMASSERT(thickness_input); + Input* thicknessobs_input=inputs->GetInput(ThicknessObsEnum);ISSMASSERT(thicknessobs_input); + Input* weights_input =inputs->GetInput(WeightsEnum);ISSMASSERT(weights_input); /* Start looping on the number of gaussian points: */ - for (ig=0; igbegin();igend();ig++){ + + gauss->GaussPoint(ig); /* Get Jacobian determinant: */ - GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss_l1l2l3); + GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss); /*Get parameters at gauss point*/ - thickness_input->GetParameterValue(&thickness, &gauss_l1l2l3[0]); - thicknessobs_input->GetParameterValue(&thicknessobs, &gauss_l1l2l3[0]); - weights_input->GetParameterValue(&weight, &gauss_l1l2l3[0]); + thickness_input->GetParameterValue(&thickness,gauss); + thicknessobs_input->GetParameterValue(&thicknessobs,gauss); + weights_input->GetParameterValue(&weight,gauss); /*compute ThicknessAbsMisfit*/ - Jelem+=0.5*pow(thickness-thicknessobs,2.0)*weight*Jdet*gauss_weight; - } - - xfree((void**)&first_gauss_area_coord); - xfree((void**)&second_gauss_area_coord); - xfree((void**)&third_gauss_area_coord); - xfree((void**)&gauss_weights); - - /*Return: */ + Jelem+=0.5*pow(thickness-thicknessobs,2.0)*weight*Jdet*gauss->weight; + } + + /* clean up and Return: */ + delete gauss; return Jelem; } @@ -1885,11 +1828,6 @@ /* gaussian points: */ - int num_gauss,ig; - double* first_gauss_area_coord = NULL; - double* second_gauss_area_coord = NULL; - double* third_gauss_area_coord = NULL; - double* gauss_weights = NULL; - double gauss_weight; - double gauss_l1l2l3[3]; + int ig; + GaussTria *gauss=NULL; /* parameters: */ @@ -1953,31 +1891,17 @@ } - /* Get gaussian points and weights (make this a statically initialized list of points? fstd): */ - GaussLegendreTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2); - /* Start looping on the number of gaussian points: */ - for (ig=0; igbegin();igend(); ig++){ + + gauss->GaussPoint(ig); + + GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss); + TriaRef::GetParameterValue(&misfit, &misfit_list[0],gauss); + Jelem+=misfit*Jdet*gauss->weight; + } + + /*clean up and Return: */ + delete gauss; return Jelem; } Index: /issm/trunk/src/c/objects/Elements/Tria.h =================================================================== --- /issm/trunk/src/c/objects/Elements/Tria.h (revision 5634) +++ /issm/trunk/src/c/objects/Elements/Tria.h (revision 5635) @@ -75,5 +75,4 @@ void CreateKMatrix(Mat Kgg); void CreatePVector(Vec pg); - void* GetMatPar(); void GetNodes(void** nodes); bool GetOnBed();