Index: /issm/trunk/src/c/objects/Elements/Penta.cpp =================================================================== --- /issm/trunk/src/c/objects/Elements/Penta.cpp (revision 5920) +++ /issm/trunk/src/c/objects/Elements/Penta.cpp (revision 5921) @@ -2984,25 +2984,15 @@ void Penta::CreatePVectorAdjointMacAyeal(Vec p_g){ - int i; - Tria* tria=NULL; - - /*If on water, skip: */ - if(IsOnWater()) return; - - /*Bail out if this element if: - * -> Non MacAyeal and not on the surface - * -> MacAyeal(2d model) and not on bed) */ - if (!IsOnBed()){ - return; - } - else{ - - /*This element should be collapsed into a tria element at its base. Create this tria element, - * and compute pe_g*/ - tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria (lower face). - tria->CreatePVectorAdjointHoriz(p_g); - delete tria->matice; delete tria; - return; - } + if (!IsOnBed() || IsOnWater()) return; + + /*Call Tria function*/ + Tria* tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria. + ElementVector* pe=tria->CreatePVectorAdjointHoriz(); + delete tria->matice; delete tria; + pe->AddToGlobal(p_g,NULL); + delete pe; + + /*clean up and return*/ + return; } /*}}}*/ @@ -3010,23 +3000,15 @@ void Penta::CreatePVectorAdjointPattyn(Vec p_g){ - int i; - Tria* tria=NULL; - - /*If on water, skip: */ - if(IsOnWater()) return; - - /*Bail out if this element if: - * -> Non MacAyeal and not on the surface - * -> MacAyeal(2d model) and not on bed) */ - if (!IsOnSurface()){ - return; - } - else{ - - tria=(Tria*)SpawnTria(3,4,5); //grids 3, 4 and 5 make the new tria (upper face). - tria->CreatePVectorAdjointHoriz(p_g); - delete tria->matice; delete tria; - return; - } + if (!IsOnSurface() || IsOnWater()) return; + + /*Call Tria function*/ + Tria* tria=(Tria*)SpawnTria(3,4,5); //grids 3, 4 and 5 make the new tria (upper face). + ElementVector* pe=tria->CreatePVectorAdjointHoriz(); + delete tria->matice; delete tria; + pe->AddToGlobal(p_g,NULL); + delete pe; + + /*clean up and return*/ + return; } /*}}}*/ @@ -3371,28 +3353,18 @@ void Penta::CreatePVectorDiagnosticMacAyeal(Vec pg){ - /*Spawning: */ - Tria* tria=NULL; - - /*If on water, skip load: */ - if(IsOnWater())return; - - /*Figure out if this pentaelem is Macayeal. If so, then bailout, except if it is at the + /*Figure out if this penta is collapsed. If so, then bailout, except if it is at the bedrock, in which case we spawn a tria element using the 3 first grids, and use it to build - the load vector. */ - - if (!IsOnBed()){ - /*This element should be collapsed, but this element is not on the bedrock, therefore all its - * dofs have already been frozen! Do nothing: */ - return; - } - else{ - - /*This element should be collapsed into a tria element at its base. Create this tria element, - *and use its CreatePVector functionality to return an elementary load vector: */ - tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria. - tria->CreatePVector(pg,NULL); - delete tria->matice; delete tria; - return; - } + the stiffness matrix. */ + if (!IsOnBed() || IsOnWater()) return; + + /*Call Tria function*/ + Tria* tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria. + ElementVector* pe=tria->CreatePVectorDiagnosticMacAyeal(); + delete tria->matice; delete tria; + pe->AddToGlobal(pg,NULL); + delete pe; + + /*clean up and return*/ + return; } /*}}}*/ Index: /issm/trunk/src/c/objects/Elements/Tria.cpp =================================================================== --- /issm/trunk/src/c/objects/Elements/Tria.cpp (revision 5920) +++ /issm/trunk/src/c/objects/Elements/Tria.cpp (revision 5921) @@ -751,8 +751,12 @@ switch(analysis_type){ case DiagnosticHorizAnalysisEnum: - CreatePVectorDiagnosticMacAyeal( pg,pf); + pe=CreatePVectorDiagnosticMacAyeal(); + pe->AddToGlobal(pg,pf); + delete pe; break; case AdjointHorizAnalysisEnum: - CreatePVectorAdjointHoriz( pg); + pe=CreatePVectorAdjointHoriz(); + pe->AddToGlobal(pg,pf); + delete pe; break; case DiagnosticHutterAnalysisEnum: @@ -3666,5 +3670,5 @@ /*}}}*/ /*FUNCTION Tria::CreatePVectorDiagnosticMacAyeal {{{1*/ -void Tria::CreatePVectorDiagnosticMacAyeal( Vec pg, Vec pf){ +ElementVector* Tria::CreatePVectorDiagnosticMacAyeal(){ /*Constants*/ @@ -3678,20 +3682,17 @@ double slope[2]; double l1l2l3[3]; - double pe_g[numdof]={0.0}; double pe_g_gaussian[numdof]; GaussTria* gauss=NULL; - ElementVector* pe=NULL; - - /*First, if we are on water, return empty vector: */ - if(IsOnWater()) return; - - /*Retrieve all Inputs and parameters: */ + + /*Initialize Element vector and return if necessary*/ + if(IsOnWater()) return NULL; + ElementVector* pe=NewElementVector(nodes,NUMVERTICES,this->parameters); + + /*Retrieve all inputs and parameters*/ + GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES); inputs->GetParameterValue(&drag_type,DragTypeEnum); 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); - - /*Get node coordinates and dof list: */ - GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES); /* Start looping on the number of gaussian points: */ @@ -3728,17 +3729,10 @@ } } - - for( i=0; iparameters,MacAyealApproximationEnum); - pe->AddValues(&pe_g[0]); - pe->AddToGlobal(pg,pf); - - /*Free ressources:*/ - delete pe; + for(i=0;ivalues[i]+=pe_g_gaussian[i]; + } + + /*Clean up and return*/ delete gauss; - - + return pe; } /*}}}*/ @@ -3798,14 +3792,9 @@ /*}}}*/ /*FUNCTION Tria::CreatePVectorAdjointHoriz{{{1*/ -void Tria::CreatePVectorAdjointHoriz(Vec p_g){ - - int i; - - /* node data: */ - const int numdof=2*NUMVERTICES; +ElementVector* Tria::CreatePVectorAdjointHoriz(void){ + + /*Constants*/ + const int numdof=NDOF2*NUMVERTICES; double xyz_list[NUMVERTICES][3]; - int* doflist=NULL; - - /* grid data: */ double vx_list[NUMVERTICES]; double vy_list[NUMVERTICES]; @@ -3815,25 +3804,13 @@ double duy_list[NUMVERTICES]; double weights_list[NUMVERTICES]; - - /* gaussian points: */ + int i; int ig; GaussTria* gauss=NULL; - - /* parameters: */ double obs_velocity_mag,velocity_mag; double dux,duy; double meanvel, epsvel; - - /*element vector : */ - double pe_g[numdof]={0.0}; double pe_g_gaussian[numdof]; - - /* Jacobian: */ double Jdet; - - /*nodal functions: */ double l1l2l3[3]; - - /*relative and algorithmic fitting: */ double scalex=0; double scaley=0; @@ -3842,13 +3819,12 @@ int response; - /*retrieve some parameters: */ + /*Initialize Element vector and return if necessary*/ + if(IsOnWater()) return NULL; + ElementVector* pe=NewElementVector(nodes,NUMVERTICES,this->parameters); + + /*Retrieve all inputs and parameters*/ + GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES); this->parameters->FindParam(&meanvel,MeanVelEnum); this->parameters->FindParam(&epsvel,EpsVelEnum); - - /* Get node coordinates and dof list: */ - GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES); - GetDofList(&doflist,NoneApproximationEnum,GsetEnum); - - /* Recover input data: */ GetParameterListOnVertices(&obs_vx_list[0],VxObsEnum); GetParameterListOnVertices(&obs_vy_list[0],VyObsEnum); @@ -3856,5 +3832,4 @@ GetParameterListOnVertices(&vy_list[0],VyEnum); GetParameterListOnVertices(&weights_list[0],WeightsEnum); - inputs->GetParameterValue(&response,CmResponseEnum); if(response==SurfaceAverageVelMisfitEnum){ @@ -3986,17 +3961,10 @@ gauss->GaussPoint(ig); - /* Get Jacobian determinant: */ GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss); - - /* Get nodal functions value at gaussian point:*/ 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); TriaRef::GetParameterValue(&duy, &duy_list[0],gauss); - /*compute Du*/ for (i=0;iweight*l1l2l3[i]; @@ -4004,16 +3972,10 @@ } - /*Add pe_g_gaussian vector to pe_g: */ - for( i=0; ivalues[i]+=pe_g_gaussian[i]; + } /*Clean up and return*/ delete gauss; - xfree((void**)&doflist); + return pe; } /*}}}*/ Index: /issm/trunk/src/c/objects/Elements/Tria.h =================================================================== --- /issm/trunk/src/c/objects/Elements/Tria.h (revision 5920) +++ /issm/trunk/src/c/objects/Elements/Tria.h (revision 5921) @@ -143,6 +143,6 @@ void CreatePVectorBalancedvelocities(Vec pg); void CreatePVectorDiagnosticBaseVert(Vec pg); - void CreatePVectorDiagnosticMacAyeal(Vec pg,Vec pf); - void CreatePVectorAdjointHoriz(Vec pg); + ElementVector* CreatePVectorDiagnosticMacAyeal(void); + ElementVector* CreatePVectorAdjointHoriz(void); void CreatePVectorAdjointStokes(Vec pg); void CreatePVectorAdjointBalancedthickness(Vec pg);