Index: /issm/trunk/src/c/objects/Elements/Tria.cpp =================================================================== --- /issm/trunk/src/c/objects/Elements/Tria.cpp (revision 5638) +++ /issm/trunk/src/c/objects/Elements/Tria.cpp (revision 5639) @@ -4677,11 +4677,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: */ @@ -4712,28 +4707,23 @@ 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); - /* 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 nodal functions value at gaussian point:*/ - GetNodalFunctions(l1l2l3, gauss_l1l2l3); + GetNodalFunctions(l1l2l3, 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); for (i=0;iweight*l1l2l3[i]; } @@ -4747,9 +4737,6 @@ VecSetValues(p_g,numdof,doflist,(const double*)pe_g,ADD_VALUES); - /*Free ressources:*/ - xfree((void**)&first_gauss_area_coord); - xfree((void**)&second_gauss_area_coord); - xfree((void**)&third_gauss_area_coord); - xfree((void**)&gauss_weights); + /*Clean up and return*/ + delete gauss; xfree((void**)&doflist); } @@ -4777,11 +4764,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: */ @@ -4945,31 +4927,26 @@ } - /* 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); /* Get Jacobian determinant: */ - GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss_l1l2l3); + GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss); /* Get nodal functions value at gaussian point:*/ - GetNodalFunctions(l1l2l3, gauss_l1l2l3); + GetNodalFunctions(l1l2l3, gauss); /*Build due_g_gaussian vector: we have three cases here, according to which type of misfit we are using. */ /*Compute absolute(x/y) at gaussian point: */ - TriaRef::GetParameterValue(&dux, &dux_list[0],gauss_l1l2l3); - TriaRef::GetParameterValue(&duy, &duy_list[0],gauss_l1l2l3); + TriaRef::GetParameterValue(&dux, &dux_list[0],gauss); + TriaRef::GetParameterValue(&duy, &duy_list[0],gauss); /*compute Du*/ for (i=0;iweight*l1l2l3[i]; + pe_g_gaussian[i*NDOF2+1]=duy*Jdet*gauss->weight*l1l2l3[i]; } @@ -4983,9 +4960,6 @@ VecSetValues(p_g,numdof,doflist,(const double*)pe_g,ADD_VALUES); - /*Free ressources:*/ - xfree((void**)&first_gauss_area_coord); - xfree((void**)&second_gauss_area_coord); - xfree((void**)&third_gauss_area_coord); - xfree((void**)&gauss_weights); + /*Clean up and return*/ + delete gauss; xfree((void**)&doflist); } @@ -5013,11 +4987,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: */ @@ -5181,31 +5150,26 @@ } - /* 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); /* Get Jacobian determinant: */ - GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss_l1l2l3); + GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss); /* Get nodal functions value at gaussian point:*/ - GetNodalFunctions(l1l2l3, gauss_l1l2l3); + GetNodalFunctions(l1l2l3, gauss); /*Build due_g_gaussian vector: we have three cases here, according to which type of misfit we are using. */ /*Compute absolute(x/y) at gaussian point: */ - TriaRef::GetParameterValue(&dux, &dux_list[0],gauss_l1l2l3); - TriaRef::GetParameterValue(&duy, &duy_list[0],gauss_l1l2l3); + TriaRef::GetParameterValue(&dux, &dux_list[0],gauss); + TriaRef::GetParameterValue(&duy, &duy_list[0],gauss); /*compute Du*/ for (i=0;iweight*l1l2l3[i]; + pe_g[i*NDOF4+1]+=duy*Jdet*gauss->weight*l1l2l3[i]; } } @@ -5214,9 +5178,6 @@ VecSetValues(p_g,numdof,doflist,(const double*)pe_g,ADD_VALUES); - /*Free ressources:*/ - xfree((void**)&first_gauss_area_coord); - xfree((void**)&second_gauss_area_coord); - xfree((void**)&third_gauss_area_coord); - xfree((void**)&gauss_weights); + /*Clean up and return*/ + delete gauss; xfree((void**)&doflist); } @@ -5244,5 +5205,5 @@ Input* slopey_input=NULL; Input* thickness_input=NULL; - double gauss[numvertices][numvertices] = {{1,0,0},{0,1,0},{0,0,1}}; //TO BE DELETED + GaussTria* gauss=NULL; /*recover some inputs: */ @@ -5266,10 +5227,14 @@ /*Spawn 3 sing elements: */ - for(i=0;i<3;i++){ + gauss=new GaussTria(); + for(i=0;iGaussVertex(i); + connectivity=nodes[i]->GetConnectivity(); - thickness_input->GetParameterValue(&thickness,&gauss[i][0]); - slopex_input->GetParameterValue(&slope[0],&gauss[i][0]); - slopey_input->GetParameterValue(&slope[1],&gauss[i][0]); + thickness_input->GetParameterValue(&thickness,gauss); + slopex_input->GetParameterValue(&slope[0],gauss); + slopey_input->GetParameterValue(&slope[1],gauss); slope2=pow(slope[0],2)+pow(slope[1],2); @@ -5287,6 +5252,7 @@ VecSetValues(pg,numdofs,doflist,(const double*)pe_g,ADD_VALUES); - - /*Free ressources:*/ + + /*Clean up and return*/ + delete gauss; xfree((void**)&doflist); } @@ -5294,5 +5260,4 @@ /*FUNCTION Tria::CreatePVectorPrognostic_CG {{{1*/ void Tria::CreatePVectorPrognostic_CG(Vec pg){ - /* local declarations */ @@ -5308,11 +5273,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; /* matrix */ @@ -5339,7 +5299,4 @@ GetDofList(&doflist); - /* 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); - /*Retrieve all inputs we will be needing: */ accumulation_input=inputs->GetInput(AccumulationRateEnum); @@ -5348,24 +5305,22 @@ /* Start looping on the number of gaussian points: */ - for (ig=0; igbegin();igend();ig++){ + + gauss->GaussPoint(ig); /* Get Jacobian determinant: */ - GetJacobianDeterminant2d(&Jdettria, &xyz_list[0][0],gauss_l1l2l3); + GetJacobianDeterminant2d(&Jdettria, &xyz_list[0][0],gauss); /*Get L matrix: */ - GetL(&L[0], &xyz_list[0][0], gauss_l1l2l3,numberofdofspernode); + GetL(&L[0], &xyz_list[0][0], gauss,numberofdofspernode); /* Get accumulation, melting and thickness at gauss point */ - accumulation_input->GetParameterValue(&accumulation_g, &gauss_l1l2l3[0]); - melting_input->GetParameterValue(&melting_g, &gauss_l1l2l3[0]); - thickness_input->GetParameterValue(&thickness_g, &gauss_l1l2l3[0]); + accumulation_input->GetParameterValue(&accumulation_g,gauss); + melting_input->GetParameterValue(&melting_g,gauss); + thickness_input->GetParameterValue(&thickness_g,gauss); /* Add value into pe_g: */ - for( i=0; iweight*(thickness_g+dt*(accumulation_g-melting_g))*L[i]; } // for (ig=0; igGetInput(AccumulationRateEnum); @@ -5436,24 +5380,22 @@ /* Start looping on the number of gaussian points: */ - for (ig=0; igbegin();igend();ig++){ + + gauss->GaussPoint(ig); /* Get Jacobian determinant: */ - GetJacobianDeterminant2d(&Jdettria, &xyz_list[0][0],gauss_l1l2l3); + GetJacobianDeterminant2d(&Jdettria, &xyz_list[0][0],gauss); /*Get L matrix: */ - GetL(&L[0], &xyz_list[0][0], gauss_l1l2l3,numberofdofspernode); + GetL(&L[0], &xyz_list[0][0], gauss,numberofdofspernode); /* Get accumulation, melting and thickness at gauss point */ - accumulation_input->GetParameterValue(&accumulation_g, &gauss_l1l2l3[0]); - melting_input->GetParameterValue(&melting_g, &gauss_l1l2l3[0]); - thickness_input->GetParameterValue(&thickness_g, &gauss_l1l2l3[0]); + accumulation_input->GetParameterValue(&accumulation_g,gauss); + melting_input->GetParameterValue(&melting_g,gauss); + thickness_input->GetParameterValue(&thickness_g,gauss); /* Add value into pe_g: */ - for( i=0; iweight*(thickness_g+dt*(accumulation_g-melting_g))*L[i]; } // for (ig=0; igGetParameterDerivativeValue(&slope[0],&xyz_list[0][0],&gauss_l1l2l3[0]); + gauss=new GaussTria(2); + for(ig=gauss->begin();igend();ig++){ + + gauss->GaussPoint(ig); + + slope_input->GetParameterDerivativeValue(&slope[0],&xyz_list[0][0],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); /*Build pe_g_gaussian vector: */ if ( (analysis_type==SurfaceSlopeXAnalysisEnum) || (analysis_type==BedSlopeXAnalysisEnum)){ - for(i=0;iweight*slope[0]*l1l2l3[i]; } if ( (analysis_type==SurfaceSlopeYAnalysisEnum) || (analysis_type==BedSlopeYAnalysisEnum)){ - for(i=0;iweight*slope[1]*l1l2l3[i]; } @@ -5556,9 +5485,6 @@ VecSetValues(pg,numdof,doflist,(const double*)pe_g,ADD_VALUES); - /*Free ressources:*/ - xfree((void**)&first_gauss_area_coord); - xfree((void**)&second_gauss_area_coord); - xfree((void**)&third_gauss_area_coord); - xfree((void**)&gauss_weights); + /*Clean up and return*/ + delete gauss; xfree((void**)&doflist); } @@ -5588,11 +5514,6 @@ /* gaussian points: */ - int num_area_gauss,ig; - double* gauss_weights = NULL; - double* first_gauss_area_coord = NULL; - double* second_gauss_area_coord = NULL; - double* third_gauss_area_coord = NULL; - double gauss_weight; - double gauss_coord[3]; + int ig; + GaussTria* gauss=NULL; /*matrices: */ @@ -5625,28 +5546,25 @@ /* Ice/ocean heat exchange flux on ice shelf base */ - GaussLegendreTria (&num_area_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2); - /*Retrieve all inputs we will be needing: */ pressure_input=inputs->GetInput(PressureEnum); /* Start looping on the number of gauss 2d (nodes on the bedrock) */ - for (ig=0; igbegin();igend();ig++){ + + gauss->GaussPoint(ig); //Get the Jacobian determinant - GetJacobianDeterminant3d(&Jdet, &xyz_list[0][0], gauss_coord); + GetJacobianDeterminant3d(&Jdet, &xyz_list[0][0], gauss); /*Get nodal functions values: */ - GetNodalFunctions(&l1l2l3[0], gauss_coord); + GetNodalFunctions(&l1l2l3[0], gauss); /*Get geothermal flux and basal friction */ - pressure_input->GetParameterValue(&pressure, &gauss_coord[0]); + pressure_input->GetParameterValue(&pressure,gauss); t_pmp=meltingpoint-beta*pressure; /*Calculate scalar parameter*/ - scalar_ocean=gauss_weight*Jdet*rho_water*mixed_layer_capacity*thermal_exchange_velocity*(t_pmp)/(heatcapacity*rho_ice); + scalar_ocean=gauss->weight*Jdet*rho_water*mixed_layer_capacity*thermal_exchange_velocity*(t_pmp)/(heatcapacity*rho_ice); if(dt){ scalar_ocean=dt*scalar_ocean; @@ -5661,9 +5579,6 @@ VecSetValues(pg,numdof,doflist,(const double*)P_terms,ADD_VALUES); - /*Free ressources:*/ - xfree((void**)&first_gauss_area_coord); - xfree((void**)&second_gauss_area_coord); - xfree((void**)&third_gauss_area_coord); - xfree((void**)&gauss_weights); + /*Clean up and return*/ + delete gauss; xfree((void**)&doflist); } @@ -5693,5 +5608,5 @@ double geothermalflux_value; double alpha2_list[numvertices]; //TO BE DELETED - double gauss[numvertices][numvertices] = {{1,0,0},{0,1,0},{0,0,1}}; //TO BE DELETED + double gaussgrids[numvertices][numvertices] = {{1,0,0},{0,1,0},{0,0,1}}; //TO BE DELETED double vx_list[numvertices]; //TO BE DELETED double vy_list[numvertices]; //TO BE DELETED @@ -5699,11 +5614,6 @@ /* gaussian points: */ - int num_area_gauss,ig; - double* gauss_weights = NULL; - double* first_gauss_area_coord = NULL; - double* second_gauss_area_coord = NULL; - double* third_gauss_area_coord = NULL; - double gauss_weight; - double gauss_coord[3]; + int ig; + GaussTria* gauss=NULL; /*matrices: */ @@ -5748,8 +5658,8 @@ /*COMPUT alpha2_list and basalfriction_list (TO BE DELETED)*/ for(i=0;iGetAlpha2(&alpha2_list[i],&gauss[i][0],VxEnum,VyEnum,VzEnum); //TO BE DELETED - } - vx_input->GetParameterValues(&vx_list[0],&gauss[0][0],3); //TO BE DELETED - vy_input->GetParameterValues(&vy_list[0],&gauss[0][0],3); //TO BE DELETED + friction->GetAlpha2(&alpha2_list[i],&gaussgrids[i][0],VxEnum,VyEnum,VzEnum); //TO BE DELETED + } + vx_input->GetParameterValues(&vx_list[0],&gaussgrids[0][0],3); //TO BE DELETED + vy_input->GetParameterValues(&vy_list[0],&gaussgrids[0][0],3); //TO BE DELETED for(i=0;ibegin();igend();ig++){ + + gauss->GaussPoint(ig); //Get the Jacobian determinant - GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0], gauss_coord); + GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0], gauss); /*Get nodal functions values: */ - GetNodalFunctions(&l1l2l3[0], gauss_coord); + GetNodalFunctions(&l1l2l3[0], gauss); /*Get geothermal flux and basal friction */ - geothermalflux_input->GetParameterValue(&geothermalflux_value, &gauss_coord[0]); + geothermalflux_input->GetParameterValue(&geothermalflux_value,gauss); /*Friction: */ //friction->GetAlpha2(&alpha2,&gauss_coord[0],VxEnum,VyEnum,VzEnum); - TriaRef::GetParameterValue(&basalfriction,&basalfriction_list[0],gauss_coord); // TO BE DELETED + TriaRef::GetParameterValue(&basalfriction,&basalfriction_list[0],gauss); // TO BE DELETED /*Calculate scalar parameter*/ - scalar=gauss_weight*Jdet*(basalfriction+geothermalflux_value)/(heatcapacity*rho_ice); + scalar=gauss->weight*Jdet*(basalfriction+geothermalflux_value)/(heatcapacity*rho_ice); if(dt){ scalar=dt*scalar; @@ -5794,11 +5702,8 @@ VecSetValues(pg,numdof,doflist,(const double*)P_terms,ADD_VALUES); - /*Free ressources:*/ - xfree((void**)&first_gauss_area_coord); - xfree((void**)&second_gauss_area_coord); - xfree((void**)&third_gauss_area_coord); - xfree((void**)&gauss_weights); + /*Clean up and return*/ + delete gauss; + delete friction; xfree((void**)&doflist); - delete friction; } /*}}}*/