Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp (revision 15725) +++ /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp (revision 15726) @@ -5286,12 +5286,190 @@ ElementVector* Penta::CreatePVectorAdjointHO(void){ + + /*Nothing to be done if not on surface*/ if (!IsOnSurface()) return NULL; - /*Call Tria function*/ - Tria* tria=(Tria*)SpawnTria(3,4,5); //nodes 3, 4 and 5 make the new tria (upper face). - ElementVector* pe=tria->CreatePVectorAdjointHoriz(); - delete tria->material; delete tria; - - /*clean up and return*/ + /*Intermediaries */ + int i,j,resp; + int *responses=NULL; + int num_responses; + IssmDouble Jdet2d; + IssmDouble obs_velocity_mag,velocity_mag; + IssmDouble dux,duy; + IssmDouble epsvel=2.220446049250313e-16; + IssmDouble meanvel=3.170979198376458e-05; /*1000 m/yr*/ + IssmDouble scalex=0.,scaley=0.,scale=0.,S=0.; + IssmDouble vx,vy,vxobs,vyobs,weight; + IssmDouble xyz_list[NUMVERTICES][3]; + IssmDouble xyz_list_tria[NUMVERTICES2D][3]; + + /*Fetch number of nodes and dof for this finite element*/ + int numnodes = this->NumberofNodes(); + + /*Initialize Element vector*/ + ElementVector* pe = new ElementVector(nodes,numnodes,this->parameters); + IssmDouble* basis = xNew(numnodes); + + /*Retrieve all inputs and parameters*/ + GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); + this->parameters->FindParam(&num_responses,InversionNumCostFunctionsEnum); + this->parameters->FindParam(&responses,NULL,NULL,StepResponsesEnum); + Input* weights_input=inputs->GetInput(InversionCostFunctionsCoefficientsEnum); _assert_(weights_input); + Input* vx_input =inputs->GetInput(VxEnum); _assert_(vx_input); + Input* vy_input =inputs->GetInput(VyEnum); _assert_(vy_input); + Input* vxobs_input =inputs->GetInput(InversionVxObsEnum); _assert_(vxobs_input); + Input* vyobs_input =inputs->GetInput(InversionVyObsEnum); _assert_(vyobs_input); + + for(i=0;iGetInputValue(&S,SurfaceAreaEnum); break; + } + } + + /* Start looping on the number of gaussian points: */ + GaussPenta* gauss=new GaussPenta(3,4,5,4); + for(int ig=gauss->begin();igend();ig++){ + + gauss->GaussPoint(ig); + + /* Get Jacobian determinant: */ + GetTriaJacobianDeterminant(&Jdet2d,&xyz_list_tria[0][0],gauss); + + /*Get all parameters at gaussian point*/ + vx_input->GetInputValue(&vx,gauss); + vy_input->GetInputValue(&vy,gauss); + vxobs_input->GetInputValue(&vxobs,gauss); + vyobs_input->GetInputValue(&vyobs,gauss); + GetNodalFunctions(basis, gauss); + + /*Loop over all requested responses*/ + for(resp=0;respGetInputValue(&weight,gauss,resp); + + switch(responses[resp]){ + case SurfaceAbsVelMisfitEnum: + /* + * 1 [ 2 2 ] + * J = --- | (u - u ) + (v - v ) | + * 2 [ obs obs ] + * + * dJ + * DU = - -- = (u - u ) + * du obs + */ + for(i=0;ivalues[i*NDOF2+0]+=dux*weight*Jdet2d*gauss->weight*basis[i]; + pe->values[i*NDOF2+1]+=duy*weight*Jdet2d*gauss->weight*basis[i]; + } + break; + case SurfaceRelVelMisfitEnum: + /* + * 1 [ \bar{v}^2 2 \bar{v}^2 2 ] + * J = --- | ------------- (u - u ) + ------------- (v - v ) | + * 2 [ (u + eps)^2 obs (v + eps)^2 obs ] + * obs obs + * + * dJ \bar{v}^2 + * DU = - -- = ------------- (u - u ) + * du (u + eps)^2 obs + * obs + */ + for(i=0;ivalues[i*NDOF2+0]+=dux*weight*Jdet2d*gauss->weight*basis[i]; + pe->values[i*NDOF2+1]+=duy*weight*Jdet2d*gauss->weight*basis[i]; + } + break; + case SurfaceLogVelMisfitEnum: + /* + * [ vel + eps ] 2 + * J = 4 \bar{v}^2 | log ( ----------- ) | + * [ vel + eps ] + * obs + * + * dJ 2 * log(...) + * DU = - -- = - 4 \bar{v}^2 ------------- u + * du vel^2 + eps + * + */ + for(i=0;ivalues[i*NDOF2+0]+=dux*weight*Jdet2d*gauss->weight*basis[i]; + pe->values[i*NDOF2+1]+=duy*weight*Jdet2d*gauss->weight*basis[i]; + } + break; + case SurfaceAverageVelMisfitEnum: + /* + * 1 2 2 + * J = --- sqrt( (u - u ) + (v - v ) ) + * S obs obs + * + * dJ 1 1 + * DU = - -- = - --- ----------- * 2 (u - u ) + * du S 2 sqrt(...) obs + */ + for(i=0;ivalues[i*NDOF2+0]+=dux*weight*Jdet2d*gauss->weight*basis[i]; + pe->values[i*NDOF2+1]+=duy*weight*Jdet2d*gauss->weight*basis[i]; + } + break; + case SurfaceLogVxVyMisfitEnum: + /* + * 1 [ |u| + eps 2 |v| + eps 2 ] + * J = --- \bar{v}^2 | log ( ----------- ) + log ( ----------- ) | + * 2 [ |u |+ eps |v |+ eps ] + * obs obs + * dJ 1 u 1 + * DU = - -- = - \bar{v}^2 log(u...) --------- ---- ~ - \bar{v}^2 log(u...) ------ + * du |u| + eps |u| u + eps + */ + for(i=0;ivalues[i*NDOF2+0]+=dux*weight*Jdet2d*gauss->weight*basis[i]; + pe->values[i*NDOF2+1]+=duy*weight*Jdet2d*gauss->weight*basis[i]; + } + break; + case DragCoefficientAbsGradientEnum: + /*Nothing in P vector*/ + break; + case ThicknessAbsGradientEnum: + /*Nothing in P vector*/ + break; + case ThicknessAlongGradientEnum: + /*Nothing in P vector*/ + break; + case ThicknessAcrossGradientEnum: + /*Nothing in P vector*/ + break; + case RheologyBbarAbsGradientEnum: + /*Nothing in P vector*/ + break; + default: + _error_("response " << EnumToStringx(responses[resp]) << " not supported yet"); + } + } + } + + /*Clean up and return*/ + xDelete(basis); + xDelete(responses); + delete gauss; return pe; } @@ -5583,14 +5761,16 @@ break; case DragCoefficientAbsGradientEnum: - if (!IsOnBed()) return; - tria=(Tria*)SpawnTria(0,1,2); //nodes 0, 1 and 2 make the new tria. - tria->GradjDragGradient(gradient,resp,control_index); - delete tria->material; delete tria; + if(IsOnBed()){ + tria=(Tria*)SpawnTria(0,1,2); //nodes 0, 1 and 2 make the new tria. + tria->GradjDragGradient(gradient,resp,control_index); + delete tria->material; delete tria; + } break; case RheologyBbarAbsGradientEnum: - if (!IsOnBed()) return; - tria=(Tria*)SpawnTria(0,1,2); //nodes 0, 1 and 2 make the new tria. - tria->GradjBGradient(gradient,resp,control_index); - delete tria->material; delete tria; + if(IsOnBed()){ + tria=(Tria*)SpawnTria(0,1,2); //nodes 0, 1 and 2 make the new tria. + tria->GradjBGradient(gradient,resp,control_index); + delete tria->material; delete tria; + } break; default: Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 15725) +++ /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 15726) @@ -6373,6 +6373,6 @@ D[0][0]=D_scalar*dvxdx; D[0][1]=0.; + D[1][0]=0.; D[1][1]=D_scalar*dvydy; - D[1][0]=0.; TripleMultiply(B,2,numnodes,1, &D[0][0],2,2,0, @@ -6823,6 +6823,6 @@ D[0][0]=D_scalar*dvxdx; D[0][1]=0.; + D[1][0]=0.; D[1][1]=D_scalar*dvydy; - D[1][0]=0.; TripleMultiply(B,2,numnodes,1, &D[0][0],2,2,0, Index: /issm/trunk-jpl/src/c/classes/FemModel.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/FemModel.cpp (revision 15725) +++ /issm/trunk-jpl/src/c/classes/FemModel.cpp (revision 15726) @@ -153,5 +153,5 @@ /*Finalize PETSC for this model: */ #ifdef _HAVE_PETSC_ - _printf0_("closing Petsc\n"); + _printf0_("closing PETSc\n"); PetscFinalize(); #endif Index: /issm/trunk-jpl/src/c/modules/ModelProcessorx/DiagnosticHoriz/UpdateElementsDiagnosticHoriz.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/ModelProcessorx/DiagnosticHoriz/UpdateElementsDiagnosticHoriz.cpp (revision 15725) +++ /issm/trunk-jpl/src/c/modules/ModelProcessorx/DiagnosticHoriz/UpdateElementsDiagnosticHoriz.cpp (revision 15726) @@ -32,10 +32,11 @@ iomodel->Constant(&materials_type,MaterialsEnum); + /*Now, is the flag macayaealHO on? otherwise, do nothing: */ + if(!isSSA & !isL1L2 & !isHO & !isFS) return; + /*Fetch data needed and allocate vectors: */ iomodel->FetchData(1,FlowequationElementEquationEnum); finiteelement_list=xNewZeroInit(iomodel->numberofelements); - /*Now, is the flag macayaealHO on? otherwise, do nothing: */ - if(!isSSA & !isL1L2 & !isHO & !isFS) return; /*Do we have coupling*/