Index: /issm/trunk/src/c/objects/Elements/Tria.cpp =================================================================== --- /issm/trunk/src/c/objects/Elements/Tria.cpp (revision 5740) +++ /issm/trunk/src/c/objects/Elements/Tria.cpp (revision 5741) @@ -2605,5 +2605,5 @@ } - } // for (ig=0; igweight*(accumulation_g-melting_g-dhdt_g)*L[i]; - } // for (ig=0; igweight*(accumulation_g-melting_g-dhdt_g)*L[i]; - } // for (ig=0; igweight*(accumulation_g-melting_g)*L[i]; - } // for (ig=0; igGetParameterValue(&onwater,ElementOnWaterEnum); inputs->GetParameterValue(&drag_type,DragTypeEnum); - thickness_input=inputs->GetInput(ThicknessEnum); - surface_input=inputs->GetInput(SurfaceEnum); - drag_input=inputs->GetInput(DragCoefficientEnum); + Input* thickness_input=inputs->GetInput(ThicknessEnum); ISSMASSERT(thickness_input); + Input* surface_input=inputs->GetInput(SurfaceEnum); ISSMASSERT(surface_input); + Input* drag_input=inputs->GetInput(DragCoefficientEnum);ISSMASSERT(drag_input); /*First, if we are on water, return empty vector: */ - if(onwater)return; + if(onwater) return; /* Get node coordinates and dof list: */ @@ -4330,31 +4310,24 @@ GetDofList(&doflist,MacAyealApproximationEnum); - - /* Get gaussian points and weights: */ - GaussLegendreTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2); /*We need higher order because our load is order 2*/ - /* Start looping on the number of gaussian points: */ - for (ig=0; igGetParameterValue(&thickness, &gauss_l1l2l3[0]); - surface_input->GetParameterDerivativeValue(&slope[0],&xyz_list[0][0],&gauss_l1l2l3[0]); + gauss=new GaussTria(2); + for(ig=gauss->begin();igend();ig++){ + + gauss->GaussPoint(ig); + + thickness_input->GetParameterValue(&thickness,gauss); + surface_input->GetParameterDerivativeValue(&slope[0],&xyz_list[0][0],gauss); /*In case we have plastic basal drag, compute plastic stress at gaussian point from k1, k2 and k3 fields in the * element itself: */ if(drag_type==1){ - drag_input->GetParameterValue(&plastic_stress, &gauss_l1l2l3[0]); + drag_input->GetParameterValue(&plastic_stress,gauss); } /* Get Jacobian determinant: */ - GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss_l1l2l3); + GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss); /*Get nodal functions: */ - GetNodalFunctions(l1l2l3, gauss_l1l2l3); + GetNodalFunctions(l1l2l3, gauss); /*Compute driving stress: */ @@ -4365,5 +4338,5 @@ for (i=0;iweight*l1l2l3[i]; } } @@ -4372,5 +4345,5 @@ for (i=0;iweight*l1l2l3[i]; } } @@ -4380,5 +4353,5 @@ for( i=0; iweight*(thickness_g+dt*(accumulation_g-melting_g))*L[i]; - } // for (ig=0; igweight*(thickness_g+dt*(accumulation_g-melting_g))*L[i]; - } // for (ig=0; igGetParameterValue(&shelf,ElementOnIceShelfEnum); inputs->GetParameterValue(&drag_type,DragTypeEnum); + Input* drag_input=inputs->GetInput(DragCoefficientEnum); ISSMASSERT(drag_input); /*retrieve some parameters: */ @@ -5838,29 +5790,24 @@ friction=new Friction("2d",inputs,matpar,analysis_type); - /* 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, 4); - /* Start looping on the number of gaussian points: */ - for (ig=0; igbegin();igend();ig++){ + + gauss->GaussPoint(ig); /*Recover alpha_complement and drag: */ - if (drag_type==2) friction->GetAlphaComplement(&alpha_complement, gauss_l1l2l3,VxEnum,VyEnum); + if (drag_type==2) friction->GetAlphaComplement(&alpha_complement, gauss,VxEnum,VyEnum); else alpha_complement=0; - inputs->GetParameterValue(&drag, &gauss_l1l2l3[0],DragCoefficientEnum); + inputs->GetParameterValue(&drag,gauss,DragCoefficientEnum); /*recover lambda mu and xi: */ - inputs->GetParameterValue(&lambda, &gauss_l1l2l3[0],AdjointxEnum); - inputs->GetParameterValue(&mu, &gauss_l1l2l3[0],AdjointyEnum); - inputs->GetParameterValue(&xi, &gauss_l1l2l3[0],AdjointzEnum); + inputs->GetParameterValue(&lambda,gauss,AdjointxEnum); + inputs->GetParameterValue(&mu, gauss,AdjointyEnum); + inputs->GetParameterValue(&xi, gauss,AdjointzEnum); /*recover vx vy and vz: */ - inputs->GetParameterValue(&vx, &gauss_l1l2l3[0],VxEnum); - inputs->GetParameterValue(&vy, &gauss_l1l2l3[0],VyEnum); - inputs->GetParameterValue(&vz, &gauss_l1l2l3[0],VzEnum); + inputs->GetParameterValue(&vx, gauss,VxEnum); + inputs->GetParameterValue(&vy, gauss,VyEnum); + inputs->GetParameterValue(&vz, gauss,VzEnum); /*Get normal vector to the bed */ @@ -5872,14 +5819,14 @@ /* Get Jacobian determinant: */ - GetJacobianDeterminant3d(&Jdet, &xyz_list[0][0],gauss_l1l2l3); + GetJacobianDeterminant3d(&Jdet, &xyz_list[0][0],gauss); /* Get nodal functions value at gaussian point:*/ - GetNodalFunctions(l1l2l3, gauss_l1l2l3); + GetNodalFunctions(l1l2l3, gauss); /*Get nodal functions derivatives*/ - GetNodalFunctionsDerivatives(&dh1dh3[0][0],&xyz_list[0][0],gauss_l1l2l3); + GetNodalFunctionsDerivatives(&dh1dh3[0][0],&xyz_list[0][0],gauss); /*Get k derivative: dk/dx */ - inputs->GetParameterDerivativeValue(&dk[0],&xyz_list[0][0],&gauss_l1l2l3[0],DragCoefficientEnum); + drag_input->GetParameterDerivativeValue(&dk[0],&xyz_list[0][0],gauss); /*Build gradje_g_gaussian vector (actually -dJ/ddrag): */ @@ -5890,8 +5837,8 @@ -mu *(2*drag*alpha_complement*(vy - vz*bed_normal[1]*bed_normal[2])) -xi *(2*drag*alpha_complement*(-vx*bed_normal[0]*bed_normal[2]-vy*bed_normal[1]*bed_normal[2])) - )*Jdet*gauss_weight*l1l2l3[i]; + )*Jdet*gauss->weight*l1l2l3[i]; //Add regularization term - grade_g_gaussian[i]+= - cm_noisedmp*Jdet*gauss_weight*(dh1dh3[0][i]*dk[0]+dh1dh3[1][i]*dk[1]); + grade_g_gaussian[i]+= - cm_noisedmp*Jdet*gauss->weight*(dh1dh3[0][i]*dk[0]+dh1dh3[1][i]*dk[1]); } @@ -5903,12 +5850,6 @@ VecSetValues(gradient,NUMVERTICES,doflist1,(const double*)grade_g,ADD_VALUES); - /*Add grade_g to the inputs of this element: */ - this->inputs->AddInput(new TriaVertexInput(GradientEnum,&grade_g[0])); - - xfree((void**)&first_gauss_area_coord); - xfree((void**)&second_gauss_area_coord); - xfree((void**)&third_gauss_area_coord); - xfree((void**)&gauss_weights); delete friction; + delete gauss; } Index: /issm/trunk/src/c/objects/Loads/Numericalflux.cpp =================================================================== --- /issm/trunk/src/c/objects/Loads/Numericalflux.cpp (revision 5740) +++ /issm/trunk/src/c/objects/Loads/Numericalflux.cpp (revision 5741) @@ -652,19 +652,20 @@ int count=0; int type; - int numberofnodes=2; - - /*output: */ + int numberofnodes; int* doflist=NULL; - /*recover type: */ - inputs->GetParameterValue(&type,TypeEnum); - /*Some checks for debugging*/ ISSMASSERT(nodes); /*How many nodes? :*/ - if(type==InternalEnum)numberofnodes=4; - else if(type==BoundaryEnum)numberofnodes=2; - else ISSMERROR("type %s not supported yet",type); + inputs->GetParameterValue(&type,TypeEnum); + switch(type){ + case InternalEnum: + numberofnodes=NUMVERTICES_INTERNAL; break; + case BoundaryEnum: + numberofnodes=NUMVERTICES_BOUNDARY; break; + default: + ISSMERROR("type %s not supported yet",type); + } /*Figure out size of doflist: */ @@ -691,5 +692,4 @@ /*Build unit outward pointing vector*/ - const int numgrids=4; double vector[2]; double norm;