Ice Sheet System Model  4.18
Code documentation
Public Member Functions
AdjointHorizAnalysis Class Reference

#include <AdjointHorizAnalysis.h>

Inheritance diagram for AdjointHorizAnalysis:
Analysis

Public Member Functions

void CreateConstraints (Constraints *constraints, IoModel *iomodel)
 
void CreateLoads (Loads *loads, IoModel *iomodel)
 
void CreateNodes (Nodes *nodes, IoModel *iomodel, bool isamr=false)
 
int DofsPerNode (int **doflist, int domaintype, int approximation)
 
void UpdateElements (Elements *elements, Inputs2 *inputs2, IoModel *iomodel, int analysis_counter, int analysis_type)
 
void UpdateParameters (Parameters *parameters, IoModel *iomodel, int solution_enum, int analysis_enum)
 
void Core (FemModel *femmodel)
 
ElementVector * CreateDVector (Element *element)
 
ElementMatrix * CreateJacobianMatrix (Element *element)
 
ElementMatrix * CreateKMatrix (Element *element)
 
ElementMatrix * CreateKMatrixFS (Element *element)
 
ElementMatrix * CreateKMatrixHO (Element *element)
 
ElementMatrix * CreateKMatrixL1L2 (Element *element)
 
ElementMatrix * CreateKMatrixSSA (Element *element)
 
ElementVector * CreatePVector (Element *element)
 
ElementVector * CreatePVectorFS (Element *element)
 
ElementVector * CreatePVectorL1L2 (Element *element)
 
ElementVector * CreatePVectorHO (Element *element)
 
ElementVector * CreatePVectorSSA (Element *element)
 
void GetSolutionFromInputs (Vector< IssmDouble > *solution, Element *element)
 
void GradientJ (Vector< IssmDouble > *gradient, Element *element, int control_type, int control_index)
 
void GradientJBbarFS (Element *element, Vector< IssmDouble > *gradient, int control_index)
 
void GradientJBbarGradient (Element *element, Vector< IssmDouble > *gradient, int control_index)
 
void GradientJBinitial (Element *element, Vector< IssmDouble > *gradient, int control_index)
 
void GradientJBbarL1L2 (Element *element, Vector< IssmDouble > *gradient, int control_index)
 
void GradientJBbarHO (Element *element, Vector< IssmDouble > *gradient, int control_index)
 
void GradientJBbarSSA (Element *element, Vector< IssmDouble > *gradient, int control_index)
 
void GradientJBFS (Element *element, Vector< IssmDouble > *gradient, int control_index)
 
void GradientJBGradient (Element *element, Vector< IssmDouble > *gradient, int control_index)
 
void GradientJBHO (Element *element, Vector< IssmDouble > *gradient, int control_index)
 
void GradientJBSSA (Element *element, Vector< IssmDouble > *gradient, int control_index)
 
void GradientJDragFS (Element *element, Vector< IssmDouble > *gradient, int control_index)
 
void GradientJDragGradient (Element *element, Vector< IssmDouble > *gradient, int control_index)
 
void GradientJDragL1L2 (Element *element, Vector< IssmDouble > *gradient, int control_index)
 
void GradientJDragHO (Element *element, Vector< IssmDouble > *gradient, int control_index)
 
void GradientJDragSSA (Element *element, Vector< IssmDouble > *gradient, int control_index)
 
void GradientJDragHydroFS (Element *element, Vector< IssmDouble > *gradient, int control_index)
 
void GradientJDragHydroL1L2 (Element *element, Vector< IssmDouble > *gradient, int control_index)
 
void GradientJDragHydroHO (Element *element, Vector< IssmDouble > *gradient, int control_index)
 
void GradientJDragHydroSSA (Element *element, Vector< IssmDouble > *gradient, int control_index)
 
void GradientJDSSA (Element *element, Vector< IssmDouble > *gradient, int control_index)
 
void InputUpdateFromSolution (IssmDouble *solution, Element *element)
 
void InputUpdateFromSolutionFS (IssmDouble *solution, Element *element)
 
void InputUpdateFromSolutionHoriz (IssmDouble *solution, Element *element)
 
void UpdateConstraints (FemModel *femmodel)
 
- Public Member Functions inherited from Analysis
virtual ~Analysis ()
 

Detailed Description

Definition at line 11 of file AdjointHorizAnalysis.h.

Member Function Documentation

◆ CreateConstraints()

void AdjointHorizAnalysis::CreateConstraints ( Constraints *  constraints,
IoModel *  iomodel 
)
virtual

Implements Analysis.

Definition at line 9 of file AdjointHorizAnalysis.cpp.

9  {/*{{{*/
10  _error_("not implemented yet");
11 }/*}}}*/

◆ CreateLoads()

void AdjointHorizAnalysis::CreateLoads ( Loads *  loads,
IoModel *  iomodel 
)
virtual

Implements Analysis.

Definition at line 12 of file AdjointHorizAnalysis.cpp.

12  {/*{{{*/
13  _error_("not implemented yet");
14 }/*}}}*/

◆ CreateNodes()

void AdjointHorizAnalysis::CreateNodes ( Nodes *  nodes,
IoModel *  iomodel,
bool  isamr = false 
)
virtual

Implements Analysis.

Definition at line 15 of file AdjointHorizAnalysis.cpp.

15  {/*{{{*/
16  _error_("not implemented yet");
17 }/*}}}*/

◆ DofsPerNode()

int AdjointHorizAnalysis::DofsPerNode ( int **  doflist,
int  domaintype,
int  approximation 
)
virtual

Implements Analysis.

Definition at line 18 of file AdjointHorizAnalysis.cpp.

18  {/*{{{*/
19  _error_("not implemented");
20 }/*}}}*/

◆ UpdateElements()

void AdjointHorizAnalysis::UpdateElements ( Elements *  elements,
Inputs2 *  inputs2,
IoModel *  iomodel,
int  analysis_counter,
int  analysis_type 
)
virtual

Implements Analysis.

Definition at line 21 of file AdjointHorizAnalysis.cpp.

21  {/*{{{*/
22  _error_("not implemented yet");
23 }/*}}}*/

◆ UpdateParameters()

void AdjointHorizAnalysis::UpdateParameters ( Parameters *  parameters,
IoModel *  iomodel,
int  solution_enum,
int  analysis_enum 
)
virtual

Implements Analysis.

Definition at line 24 of file AdjointHorizAnalysis.cpp.

24  {/*{{{*/
25  _error_("not implemented yet");
26 }/*}}}*/

◆ Core()

void AdjointHorizAnalysis::Core ( FemModel *  femmodel)
virtual

Implements Analysis.

Definition at line 29 of file AdjointHorizAnalysis.cpp.

29  {/*{{{*/
30  _error_("not implemented");
31 }/*}}}*/

◆ CreateDVector()

ElementVector * AdjointHorizAnalysis::CreateDVector ( Element *  element)
virtual

Implements Analysis.

Definition at line 32 of file AdjointHorizAnalysis.cpp.

32  {/*{{{*/
33  /*Default, return NULL*/
34  return NULL;
35 }/*}}}*/

◆ CreateJacobianMatrix()

ElementMatrix * AdjointHorizAnalysis::CreateJacobianMatrix ( Element *  element)
virtual

Implements Analysis.

Definition at line 36 of file AdjointHorizAnalysis.cpp.

36  {/*{{{*/
37 _error_("Not implemented");
38 }/*}}}*/

◆ CreateKMatrix()

ElementMatrix * AdjointHorizAnalysis::CreateKMatrix ( Element *  element)
virtual

Implements Analysis.

Definition at line 39 of file AdjointHorizAnalysis.cpp.

39  {/*{{{*/
40  int approximation;
41  element->GetInputValue(&approximation,ApproximationEnum);
42  switch(approximation){
43  case SSAApproximationEnum:
44  return CreateKMatrixSSA(element);
45  case L1L2ApproximationEnum:
46  return CreateKMatrixL1L2(element);
47  case HOApproximationEnum:
48  return CreateKMatrixHO(element);
49  case FSApproximationEnum:
50  return CreateKMatrixFS(element);
51  case NoneApproximationEnum:
52  return NULL;
53  default:
54  _error_("Approximation "<<EnumToStringx(approximation)<<" not supported");
55  }
56 }/*}}}*/

◆ CreateKMatrixFS()

ElementMatrix * AdjointHorizAnalysis::CreateKMatrixFS ( Element *  element)

Definition at line 57 of file AdjointHorizAnalysis.cpp.

57  {/*{{{*/
58 
59  /* Check if ice in element */
60  if(!element->IsIceInElement()) return NULL;
61 
62  /*Intermediaries */
63  bool incomplete_adjoint;
64  int dim,epssize;
65  IssmDouble Jdet,mu_prime;
66  IssmDouble eps1dotdphii,eps1dotdphij,eps2dotdphii,eps2dotdphij,eps3dotdphii,eps3dotdphij;
67  IssmDouble eps1[3],eps2[3],eps3[3],epsilon[5];/* epsilon=[exx,eyy,exy,exz,eyz];*/
68  IssmDouble *xyz_list = NULL;
69 
70  /*Get problem dimension*/
71  element->FindParam(&dim,DomainDimensionEnum);
72  if(dim==2) epssize = 3;
73  else epssize = 6;
74 
75  /*Fetch number of nodes and dof for this finite element*/
76  int vnumnodes = element->NumberofNodesVelocity();
77  int pnumnodes = element->NumberofNodesPressure();
78  int numdof = vnumnodes*dim + pnumnodes;
79 
80  /*Initialize Jacobian with regular FS (first part of the Gateau derivative)*/
81  element->FindParam(&incomplete_adjoint,InversionIncompleteAdjointEnum);
82  StressbalanceAnalysis* analysis = new StressbalanceAnalysis();
83  ElementMatrix* Ke=analysis->CreateKMatrix(element);
84  delete analysis;
85  if(incomplete_adjoint) return Ke;
86 
87  /*Retrieve all inputs and parameters*/
88  element->GetVerticesCoordinates(&xyz_list);
89  Input2* vx_input = element->GetInput2(VxEnum);_assert_(vx_input);
90  Input2* vy_input = element->GetInput2(VyEnum);_assert_(vy_input);
91  Input2* vz_input = NULL;
92  if(dim==3){
93  vz_input = element->GetInput2(VzEnum);
94  }
95  else{
96  _error_("Not implemented yet");
97  }
98 
99  /*Allocate dbasis*/
100  IssmDouble* dbasis = xNew<IssmDouble>(dim*vnumnodes);
101 
102  /* Start looping on the number of gaussian points: */
103  Gauss* gauss=element->NewGauss(5);
104  for(int ig=gauss->begin();ig<gauss->end();ig++){
105  gauss->GaussPoint(ig);
106 
107  element->JacobianDeterminant(&Jdet,xyz_list,gauss);
108  element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
109 
110  element->StrainRateHO(&epsilon[0],xyz_list,gauss,vx_input,vy_input);
111  element->material->ViscosityFSDerivativeEpsSquare(&mu_prime,&epsilon[0],gauss);
112  eps1[0]=epsilon[0]; eps2[0]=epsilon[2]; eps3[0]=epsilon[3];
113  eps1[1]=epsilon[2]; eps2[1]=epsilon[1]; eps3[1]=epsilon[4];
114  eps1[2]=epsilon[3]; eps2[2]=epsilon[4]; eps3[2]= -epsilon[0] -epsilon[1];
115 
116  for(int i=0;i<vnumnodes;i++){
117  for(int j=0;j<vnumnodes;j++){
118  eps1dotdphii=eps1[0]*dbasis[0*vnumnodes+i]+eps1[1]*dbasis[1*vnumnodes+i]+eps1[2]*dbasis[2*vnumnodes+i];
119  eps1dotdphij=eps1[0]*dbasis[0*vnumnodes+j]+eps1[1]*dbasis[1*vnumnodes+j]+eps1[2]*dbasis[2*vnumnodes+j];
120  eps2dotdphii=eps2[0]*dbasis[0*vnumnodes+i]+eps2[1]*dbasis[1*vnumnodes+i]+eps2[2]*dbasis[2*vnumnodes+i];
121  eps2dotdphij=eps2[0]*dbasis[0*vnumnodes+j]+eps2[1]*dbasis[1*vnumnodes+j]+eps2[2]*dbasis[2*vnumnodes+j];
122  eps3dotdphii=eps3[0]*dbasis[0*vnumnodes+i]+eps3[1]*dbasis[1*vnumnodes+i]+eps3[2]*dbasis[2*vnumnodes+i];
123  eps3dotdphij=eps3[0]*dbasis[0*vnumnodes+j]+eps3[1]*dbasis[1*vnumnodes+j]+eps3[2]*dbasis[2*vnumnodes+j];
124 
125  Ke->values[numdof*(4*i+0)+4*j+0]+=gauss->weight*Jdet*2*mu_prime*eps1dotdphij*eps1dotdphii;
126  Ke->values[numdof*(4*i+0)+4*j+1]+=gauss->weight*Jdet*2*mu_prime*eps2dotdphij*eps1dotdphii;
127  Ke->values[numdof*(4*i+0)+4*j+2]+=gauss->weight*Jdet*2*mu_prime*eps3dotdphij*eps1dotdphii;
128 
129  Ke->values[numdof*(4*i+1)+4*j+0]+=gauss->weight*Jdet*2*mu_prime*eps1dotdphij*eps2dotdphii;
130  Ke->values[numdof*(4*i+1)+4*j+1]+=gauss->weight*Jdet*2*mu_prime*eps2dotdphij*eps2dotdphii;
131  Ke->values[numdof*(4*i+1)+4*j+2]+=gauss->weight*Jdet*2*mu_prime*eps3dotdphij*eps2dotdphii;
132 
133  Ke->values[numdof*(4*i+2)+4*j+0]+=gauss->weight*Jdet*2*mu_prime*eps1dotdphij*eps3dotdphii;
134  Ke->values[numdof*(4*i+2)+4*j+1]+=gauss->weight*Jdet*2*mu_prime*eps2dotdphij*eps3dotdphii;
135  Ke->values[numdof*(4*i+2)+4*j+2]+=gauss->weight*Jdet*2*mu_prime*eps3dotdphij*eps3dotdphii;
136  }
137  }
138  }
139 
140  /*Transform Coordinate System*/
141  element->TransformStiffnessMatrixCoord(Ke,XYZEnum);
142 
143  /*Clean up and return*/
144  delete gauss;
145  xDelete<IssmDouble>(dbasis);
146  xDelete<IssmDouble>(xyz_list);
147  return Ke;
148 }/*}}}*/

◆ CreateKMatrixHO()

ElementMatrix * AdjointHorizAnalysis::CreateKMatrixHO ( Element *  element)

Definition at line 149 of file AdjointHorizAnalysis.cpp.

149  {/*{{{*/
150 
151  /* Check if ice in element */
152  if(!element->IsIceInElement()) return NULL;
153 
154  /*Intermediaries */
155  bool incomplete_adjoint;
156  IssmDouble Jdet,mu_prime;
157  IssmDouble eps1dotdphii,eps1dotdphij,eps2dotdphii,eps2dotdphij;
158  IssmDouble eps1[3],eps2[3],epsilon[5];/* epsilon=[exx,eyy,exy,exz,eyz];*/
159  IssmDouble *xyz_list = NULL;
160 
161  /*Fetch number of nodes and dof for this finite element*/
162  int numnodes = element->GetNumberOfNodes();
163 
164  /*Initialize Jacobian with regular HO (first part of the Gateau derivative)*/
165  element->FindParam(&incomplete_adjoint,InversionIncompleteAdjointEnum);
166  StressbalanceAnalysis* analysis = new StressbalanceAnalysis();
167  ElementMatrix* Ke=analysis->CreateKMatrix(element);
168  delete analysis;
169  if(incomplete_adjoint) return Ke;
170 
171  /*Retrieve all inputs and parameters*/
172  element->GetVerticesCoordinates(&xyz_list);
173  Input2* vx_input = element->GetInput2(VxEnum); _assert_(vx_input);
174  Input2* vy_input = element->GetInput2(VyEnum); _assert_(vy_input);
175 
176  /*Allocate dbasis*/
177  IssmDouble* dbasis = xNew<IssmDouble>(3*numnodes);
178 
179  /* Start looping on the number of gaussian points: */
180  Gauss* gauss=element->NewGauss(5);
181  for(int ig=gauss->begin();ig<gauss->end();ig++){
182  gauss->GaussPoint(ig);
183 
184  element->JacobianDeterminant(&Jdet,xyz_list,gauss);
185  element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
186 
187  element->StrainRateHO(&epsilon[0],xyz_list,gauss,vx_input,vy_input);
188  element->material->ViscosityHODerivativeEpsSquare(&mu_prime,&epsilon[0],gauss);
189  eps1[0]=2.*epsilon[0]+epsilon[1]; eps2[0]=epsilon[2];
190  eps1[1]=epsilon[2]; eps2[1]=epsilon[0]+2.*epsilon[1];
191  eps1[2]=epsilon[3]; eps2[2]=epsilon[4];
192 
193  for(int i=0;i<numnodes;i++){
194  for(int j=0;j<numnodes;j++){
195  eps1dotdphii=eps1[0]*dbasis[0*numnodes+i]+eps1[1]*dbasis[1*numnodes+i]+eps1[2]*dbasis[2*numnodes+i];
196  eps1dotdphij=eps1[0]*dbasis[0*numnodes+j]+eps1[1]*dbasis[1*numnodes+j]+eps1[2]*dbasis[2*numnodes+j];
197  eps2dotdphii=eps2[0]*dbasis[0*numnodes+i]+eps2[1]*dbasis[1*numnodes+i]+eps2[2]*dbasis[2*numnodes+i];
198  eps2dotdphij=eps2[0]*dbasis[0*numnodes+j]+eps2[1]*dbasis[1*numnodes+j]+eps2[2]*dbasis[2*numnodes+j];
199 
200  Ke->values[2*numnodes*(2*i+0)+2*j+0]+=gauss->weight*Jdet*2*mu_prime*eps1dotdphij*eps1dotdphii;
201  Ke->values[2*numnodes*(2*i+0)+2*j+1]+=gauss->weight*Jdet*2*mu_prime*eps2dotdphij*eps1dotdphii;
202  Ke->values[2*numnodes*(2*i+1)+2*j+0]+=gauss->weight*Jdet*2*mu_prime*eps1dotdphij*eps2dotdphii;
203  Ke->values[2*numnodes*(2*i+1)+2*j+1]+=gauss->weight*Jdet*2*mu_prime*eps2dotdphij*eps2dotdphii;
204  }
205  }
206  }
207 
208  /*Transform Coordinate System*/
209  element->TransformStiffnessMatrixCoord(Ke,XYEnum);
210 
211  /*Clean up and return*/
212  delete gauss;
213  xDelete<IssmDouble>(dbasis);
214  xDelete<IssmDouble>(xyz_list);
215  return Ke;
216 }/*}}}*/

◆ CreateKMatrixL1L2()

ElementMatrix * AdjointHorizAnalysis::CreateKMatrixL1L2 ( Element *  element)

Definition at line 217 of file AdjointHorizAnalysis.cpp.

217  {/*{{{*/
218 
219  /*Intermediaries*/
220  bool incomplete_adjoint;
221 
222  /*Initialize Jacobian with regular L1L2 (first part of the Gateau derivative)*/
223  StressbalanceAnalysis* analysis = new StressbalanceAnalysis();
224  ElementMatrix* Ke=analysis->CreateKMatrix(element);
225  delete analysis;
226 
227  /*return*/
228  element->FindParam(&incomplete_adjoint,InversionIncompleteAdjointEnum);
229  if(!incomplete_adjoint){
230  _error_("Exact adjoint not supported yet for L1L2 model");
231  }
232  return Ke;
233 }/*}}}*/

◆ CreateKMatrixSSA()

ElementMatrix * AdjointHorizAnalysis::CreateKMatrixSSA ( Element *  element)

Definition at line 234 of file AdjointHorizAnalysis.cpp.

234  {/*{{{*/
235 
236  /* Check if ice in element */
237  if(!element->IsIceInElement()) return NULL;
238 
239  /*Intermediaries*/
240  int domaintype;
241  Element* basalelement;
242 
243  /*Get basal element*/
244  element->FindParam(&domaintype,DomainTypeEnum);
245  switch(domaintype){
246  case Domain2DhorizontalEnum:
247  basalelement = element;
248  break;
249  case Domain3DEnum:
250  if(!element->IsOnBase()) return NULL;
251  basalelement = element->SpawnBasalElement();
252  break;
253  default: _error_("mesh "<<EnumToStringx(domaintype)<<" not supported yet");
254  }
255 
256  /*Intermediaries */
257  bool incomplete_adjoint;
258  IssmDouble Jdet,thickness,mu_prime;
259  IssmDouble eps1dotdphii,eps1dotdphij,eps2dotdphii,eps2dotdphij;
260  IssmDouble eps1[2],eps2[2],epsilon[3];
261  IssmDouble *xyz_list = NULL;
262 
263  /*Fetch number of nodes and dof for this finite element*/
264  int numnodes = basalelement->GetNumberOfNodes();
265 
266  /*Initialize Jacobian with regular SSA (first part of the Gateau derivative)*/
267  basalelement->FindParam(&incomplete_adjoint,InversionIncompleteAdjointEnum);
268  StressbalanceAnalysis* analysis = new StressbalanceAnalysis();
269  ElementMatrix* Ke=analysis->CreateKMatrix(element);
270  delete analysis;
271  if(incomplete_adjoint){
272  if(domaintype!=Domain2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;};
273  return Ke;
274  }
275 
276  /*Retrieve all inputs and parameters*/
277  basalelement->GetVerticesCoordinates(&xyz_list);
278  Input2* vx_input = basalelement->GetInput2(VxEnum); _assert_(vx_input);
279  Input2* vy_input = basalelement->GetInput2(VyEnum); _assert_(vy_input);
280  Input2* thickness_input = basalelement->GetInput2(ThicknessEnum); _assert_(thickness_input);
281 
282  /*Allocate dbasis*/
283  IssmDouble* dbasis = xNew<IssmDouble>(2*numnodes);
284 
285  /* Start looping on the number of gaussian points: */
286  Gauss* gauss=basalelement->NewGauss(2);
287  for(int ig=gauss->begin();ig<gauss->end();ig++){
288  gauss->GaussPoint(ig);
289 
290  basalelement->JacobianDeterminant(&Jdet,xyz_list,gauss);
291  basalelement->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
292 
293  thickness_input->GetInputValue(&thickness, gauss);
294  basalelement->StrainRateSSA(&epsilon[0],xyz_list,gauss,vx_input,vy_input);
295  basalelement->material->ViscositySSADerivativeEpsSquare(&mu_prime,&epsilon[0],gauss);
296  eps1[0]=2.*epsilon[0]+epsilon[1]; eps2[0]=epsilon[2];
297  eps1[1]=epsilon[2]; eps2[1]=epsilon[0]+2*epsilon[1];
298 
299  for(int i=0;i<numnodes;i++){
300  for(int j=0;j<numnodes;j++){
301  eps1dotdphii=eps1[0]*dbasis[0*numnodes+i]+eps1[1]*dbasis[1*numnodes+i];
302  eps1dotdphij=eps1[0]*dbasis[0*numnodes+j]+eps1[1]*dbasis[1*numnodes+j];
303  eps2dotdphii=eps2[0]*dbasis[0*numnodes+i]+eps2[1]*dbasis[1*numnodes+i];
304  eps2dotdphij=eps2[0]*dbasis[0*numnodes+j]+eps2[1]*dbasis[1*numnodes+j];
305 
306  Ke->values[2*numnodes*(2*i+0)+2*j+0]+=gauss->weight*Jdet*2*mu_prime*thickness*eps1dotdphij*eps1dotdphii;
307  Ke->values[2*numnodes*(2*i+0)+2*j+1]+=gauss->weight*Jdet*2*mu_prime*thickness*eps2dotdphij*eps1dotdphii;
308  Ke->values[2*numnodes*(2*i+1)+2*j+0]+=gauss->weight*Jdet*2*mu_prime*thickness*eps1dotdphij*eps2dotdphii;
309  Ke->values[2*numnodes*(2*i+1)+2*j+1]+=gauss->weight*Jdet*2*mu_prime*thickness*eps2dotdphij*eps2dotdphii;
310  }
311  }
312  }
313 
314  /*Transform Coordinate System*/
315  basalelement->TransformStiffnessMatrixCoord(Ke,XYEnum);
316 
317  /*Clean up and return*/
318  delete gauss;
319  xDelete<IssmDouble>(dbasis);
320  xDelete<IssmDouble>(xyz_list);
321  if(domaintype!=Domain2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;};
322  return Ke;
323 }/*}}}*/

◆ CreatePVector()

ElementVector * AdjointHorizAnalysis::CreatePVector ( Element *  element)
virtual

Implements Analysis.

Definition at line 324 of file AdjointHorizAnalysis.cpp.

324  {/*{{{*/
325 
326  int approximation;
327  element->GetInputValue(&approximation,ApproximationEnum);
328  switch(approximation){
329  case SSAApproximationEnum:
330  return CreatePVectorSSA(element);
331  case L1L2ApproximationEnum:
332  return CreatePVectorL1L2(element);
333  case HOApproximationEnum:
334  return CreatePVectorHO(element);
335  case FSApproximationEnum:
336  return CreatePVectorFS(element);
337  case NoneApproximationEnum:
338  return NULL;
339  default:
340  _error_("Approximation "<<EnumToStringx(approximation)<<" not supported");
341  }
342 }/*}}}*/

◆ CreatePVectorFS()

ElementVector * AdjointHorizAnalysis::CreatePVectorFS ( Element *  element)

Definition at line 343 of file AdjointHorizAnalysis.cpp.

343  {/*{{{*/
344 
345  /*Nothing to be done if not on surface*/
346  if(!element->IsOnSurface()) return NULL;
347 
348  /*Intermediaries */
349  int num_responses,i,domaintype;
350  IssmDouble Jdet,obs_velocity_mag,velocity_mag;
351  IssmDouble vx,vy,vxobs,vyobs,dux,duy,weight;
352  IssmDouble scalex,scaley,scale,S;
353  int *responses = NULL;
354  IssmDouble *xyz_list_top = NULL;
355 
356  /* Get domaintype*/
357  element->FindParam(&domaintype,DomainTypeEnum);
358 
359  /*Fetch number of nodes and dof for this finite element*/
360  int vnumnodes = element->NumberofNodesVelocity();
361  int pnumnodes = element->NumberofNodesPressure();
362 
363  /*Prepare coordinate system list*/
364  int* cs_list = xNew<int>(vnumnodes+pnumnodes);
365  if(domaintype==Domain2DverticalEnum) for(i=0;i<vnumnodes;i++) cs_list[i] = XYEnum;
366  else for(i=0;i<vnumnodes;i++) cs_list[i] = XYZEnum;
367  for(i=0;i<pnumnodes;i++) cs_list[vnumnodes+i] = PressureEnum;
368 
369  /*Initialize Element vector and vectors*/
370  ElementVector* pe = element->NewElementVector(FSApproximationEnum);
371  IssmDouble* vbasis = xNew<IssmDouble>(vnumnodes);
372 
373  /*Retrieve all inputs and parameters*/
374  element->GetVerticesCoordinatesTop(&xyz_list_top);
375  element->FindParam(&num_responses,InversionNumCostFunctionsEnum);
376  element->FindParam(&responses,NULL,InversionCostFunctionsEnum);
377  DatasetInput2* weights_input = element->GetDatasetInput2(InversionCostFunctionsCoefficientsEnum); _assert_(weights_input);
378  Input2* vx_input = element->GetInput2(VxEnum); _assert_(vx_input);
379  Input2* vxobs_input = element->GetInput2(InversionVxObsEnum); _assert_(vxobs_input);
380  Input2* vy_input = NULL;
381  Input2* vyobs_input = NULL;
382  if(domaintype!=Domain2DverticalEnum){
383  vy_input = element->GetInput2(VyEnum); _assert_(vy_input);
384  vyobs_input = element->GetInput2(InversionVyObsEnum); _assert_(vyobs_input);
385  }
386  IssmDouble epsvel = 2.220446049250313e-16;
387  IssmDouble meanvel = 3.170979198376458e-05; /*1000 m/yr*/
388 
389  /*Get Surface if required by one response*/
390  Input2* S_input = NULL;
391  for(int resp=0;resp<num_responses;resp++){
392  if(responses[resp]==SurfaceAverageVelMisfitEnum){
393  S_input = element->GetInput2(SurfaceAreaEnum); _assert_(S_input); break;
394  }
395  }
396 
397  /* Start looping on the number of gaussian points: */
398  Gauss* gauss=element->NewGaussTop(4);
399  for(int ig=gauss->begin();ig<gauss->end();ig++){
400  gauss->GaussPoint(ig);
401 
402  element->JacobianDeterminantTop(&Jdet,xyz_list_top,gauss);
403  element->NodalFunctionsVelocity(vbasis,gauss);
404 
405  vx_input->GetInputValue(&vx,gauss);
406  vxobs_input->GetInputValue(&vxobs,gauss);
407  if(domaintype!=Domain2DverticalEnum) {
408  vy_input->GetInputValue(&vy,gauss);
409  vyobs_input->GetInputValue(&vyobs,gauss);
410  }
411 
412  /*Loop over all requested responses*/
413  for(int resp=0;resp<num_responses;resp++){
414  weights_input->GetInputValue(&weight,gauss,responses[resp]);
415 
416  switch(responses[resp]){
417  case SurfaceAbsVelMisfitEnum:
418  /*
419  * 1 [ 2 2 ]
420  * J = --- | (u - u ) + (v - v ) |
421  * 2 [ obs obs ]
422  *
423  * dJ
424  * DU = - -- = (u - u )
425  * du obs
426  */
427  for(i=0;i<vnumnodes;i++){
428  if(domaintype!=Domain2DverticalEnum){
429  dux=vxobs-vx;
430  pe->values[i*3+0]+=dux*weight*Jdet*gauss->weight*vbasis[i];
431  duy=vyobs-vy;
432  pe->values[i*3+1]+=duy*weight*Jdet*gauss->weight*vbasis[i];
433  }
434  else{
435  dux=vxobs-vx;
436  pe->values[i*2+0]+=dux*weight*Jdet*gauss->weight*vbasis[i];
437  }
438  }
439  break;
440  case SurfaceRelVelMisfitEnum:
441  /*
442  * 1 [ \bar{v}^2 2 \bar{v}^2 2 ]
443  * J = --- | ------------- (u - u ) + ------------- (v - v ) |
444  * 2 [ (u + eps)^2 obs (v + eps)^2 obs ]
445  * obs obs
446  *
447  * dJ \bar{v}^2
448  * DU = - -- = ------------- (u - u )
449  * du (u + eps)^2 obs
450  * obs
451  */
452  for(i=0;i<vnumnodes;i++){
453  if(domaintype!=Domain2DverticalEnum){
454  scalex=pow(meanvel/(vxobs+epsvel),2); if(vxobs==0)scalex=0;
455  scaley=pow(meanvel/(vyobs+epsvel),2); if(vyobs==0)scaley=0;
456  dux=scalex*(vxobs-vx);
457  duy=scaley*(vyobs-vy);
458  pe->values[i*3+0]+=dux*weight*Jdet*gauss->weight*vbasis[i];
459  pe->values[i*3+1]+=duy*weight*Jdet*gauss->weight*vbasis[i];
460  }
461  else{
462  scalex=pow(meanvel/(vxobs+epsvel),2); if(vxobs==0)scalex=0;
463  dux=scalex*(vxobs-vx);
464  pe->values[i*2+0]+=dux*weight*Jdet*gauss->weight*vbasis[i];
465  pe->values[i*2+1]+=duy*weight*Jdet*gauss->weight*vbasis[i];
466  }
467  }
468  break;
469  case SurfaceLogVelMisfitEnum:
470  /*
471  * [ vel + eps ] 2
472  * J = 4 \bar{v}^2 | log ( ----------- ) |
473  * [ vel + eps ]
474  * obs
475  *
476  * dJ 2 * log(...)
477  * DU = - -- = - 4 \bar{v}^2 ------------- u
478  * du vel^2 + eps
479  *
480  */
481  for(i=0;i<vnumnodes;i++){
482  if(domaintype!=Domain2DverticalEnum){
483  velocity_mag =sqrt(vx*vx+vy*vy)+epsvel;
484  obs_velocity_mag=sqrt(vxobs*vxobs+vyobs*vyobs)+epsvel;
485  scale=-8.*meanvel*meanvel/(velocity_mag*velocity_mag)*log(velocity_mag/obs_velocity_mag);
486  dux=scale*vx;
487  duy=scale*vy;
488  pe->values[i*3+0]+=dux*weight*Jdet*gauss->weight*vbasis[i];
489  pe->values[i*3+1]+=duy*weight*Jdet*gauss->weight*vbasis[i];
490  }
491  else{
492  velocity_mag =fabs(vx)+epsvel;
493  obs_velocity_mag=fabs(vxobs)+epsvel;
494  scale=-8.*meanvel*meanvel/(velocity_mag*velocity_mag)*log(velocity_mag/obs_velocity_mag);
495  dux=scale*vx;
496  pe->values[i*2+0]+=dux*weight*Jdet*gauss->weight*vbasis[i];
497  }
498  }
499  break;
500  case SurfaceAverageVelMisfitEnum:
501  /*
502  * 1 2 2
503  * J = --- sqrt( (u - u ) + (v - v ) )
504  * S obs obs
505  *
506  * dJ 1 1
507  * DU = - -- = - --- ----------- * 2 (u - u )
508  * du S 2 sqrt(...) obs
509  */
510  S_input->GetInputValue(&S,gauss);
511  for(i=0;i<vnumnodes;i++){
512  if (domaintype!=Domain2DverticalEnum){
513  scale=1./(S*2*sqrt(pow(vx-vxobs,2)+pow(vy-vyobs,2))+epsvel);
514  dux=scale*(vxobs-vx);
515  duy=scale*(vyobs-vy);
516  pe->values[i*3+0]+=dux*weight*Jdet*gauss->weight*vbasis[i];
517  pe->values[i*3+1]+=duy*weight*Jdet*gauss->weight*vbasis[i];
518  }
519  else{
520  scale=1./(S*2*fabs(vx-vxobs)+epsvel);
521  dux=scale*(vxobs-vx);
522  pe->values[i*2+0]+=dux*weight*Jdet*gauss->weight*vbasis[i];
523  }
524  }
525  break;
526  case SurfaceLogVxVyMisfitEnum:
527  /*
528  * 1 [ |u| + eps 2 |v| + eps 2 ]
529  * J = --- \bar{v}^2 | log ( ----------- ) + log ( ----------- ) |
530  * 2 [ |u |+ eps |v |+ eps ]
531  * obs obs
532  * dJ 1 u 1
533  * DU = - -- = - \bar{v}^2 log(u...) --------- ---- ~ - \bar{v}^2 log(u...) ------
534  * du |u| + eps |u| u + eps
535  */
536  for(i=0;i<vnumnodes;i++){
537  if(domaintype!=Domain2DverticalEnum){
538  dux = - meanvel*meanvel * log((fabs(vx)+epsvel)/(fabs(vxobs)+epsvel)) / (vx+epsvel);
539  duy = - meanvel*meanvel * log((fabs(vy)+epsvel)/(fabs(vyobs)+epsvel)) / (vy+epsvel);
540  pe->values[i*3+0]+=dux*weight*Jdet*gauss->weight*vbasis[i];
541  pe->values[i*3+1]+=duy*weight*Jdet*gauss->weight*vbasis[i];
542  }
543  else{
544  dux = - meanvel*meanvel * log((fabs(vx)+epsvel)/(fabs(vxobs)+epsvel)) / (vx+epsvel);
545  pe->values[i*2+0]+=dux*weight*Jdet*gauss->weight*vbasis[i];
546  }
547  }
548  break;
549  case DragCoefficientAbsGradientEnum:
550  /*Nothing in P vector*/
551  break;
552  case ThicknessAbsGradientEnum:
553  /*Nothing in P vector*/
554  break;
555  case ThicknessAlongGradientEnum:
556  /*Nothing in P vector*/
557  break;
558  case ThicknessAcrossGradientEnum:
559  /*Nothing in P vector*/
560  break;
561  case RheologyBbarAbsGradientEnum:
562  /*Nothing in P vector*/
563  break;
564  case RheologyBInitialguessMisfitEnum:
565  /*Nothing in P vector*/
566  break;
567  default:
568  _error_("response " << EnumToStringx(responses[resp]) << " not supported yet");
569  }
570  }
571  }
572 
573  /*Transform coordinate system*/
574  element->TransformLoadVectorCoord(pe,cs_list);
575 
576  /*Clean up and return*/
577  xDelete<int>(cs_list);
578  xDelete<int>(responses);
579  xDelete<IssmDouble>(xyz_list_top);
580  xDelete<IssmDouble>(vbasis);
581  delete gauss;
582  return pe;
583 }/*}}}*/

◆ CreatePVectorL1L2()

ElementVector * AdjointHorizAnalysis::CreatePVectorL1L2 ( Element *  element)

Definition at line 584 of file AdjointHorizAnalysis.cpp.

584  {/*{{{*/
585 
586  /*Same as SSA*/
587  return this->CreatePVectorSSA(element);
588 }/*}}}*/

◆ CreatePVectorHO()

ElementVector * AdjointHorizAnalysis::CreatePVectorHO ( Element *  element)

Definition at line 589 of file AdjointHorizAnalysis.cpp.

589  {/*{{{*/
590 
591  /*Nothing to be done if not on surface*/
592  if(!element->IsOnSurface()) return NULL;
593 
594  /*Intermediaries */
595  int num_responses,i,domaintype;
596  IssmDouble Jdet,obs_velocity_mag,velocity_mag;
597  IssmDouble vx,vy,vxobs,vyobs,dux,duy,weight;
598  IssmDouble scalex,scaley,scale,S;
599  int *responses = NULL;
600  IssmDouble *xyz_list_top = NULL;
601 
602  /* Get domaintype*/
603  element->FindParam(&domaintype,DomainTypeEnum);
604 
605  /*Fetch number of nodes and dof for this finite element*/
606  int numnodes = element->GetNumberOfNodes();
607 
608  /*Initialize Element vector and vectors*/
609  ElementVector* pe = element->NewElementVector(HOApproximationEnum);
610  IssmDouble* basis = xNew<IssmDouble>(numnodes);
611 
612  /*Retrieve all inputs and parameters*/
613  element->GetVerticesCoordinatesTop(&xyz_list_top);
614  element->FindParam(&num_responses,InversionNumCostFunctionsEnum);
615  element->FindParam(&responses,NULL,InversionCostFunctionsEnum);
616  DatasetInput2* weights_input = element->GetDatasetInput2(InversionCostFunctionsCoefficientsEnum); _assert_(weights_input);
617  Input2* vx_input = element->GetInput2(VxEnum); _assert_(vx_input);
618  Input2* vxobs_input = element->GetInput2(InversionVxObsEnum); _assert_(vxobs_input);
619  Input2* vy_input=NULL;
620  Input2* vyobs_input=NULL;
621  if(domaintype!=Domain2DverticalEnum){
622  vy_input = element->GetInput2(VyEnum); _assert_(vy_input);
623  vyobs_input = element->GetInput2(InversionVyObsEnum); _assert_(vyobs_input);
624  }
625  IssmDouble epsvel = 2.220446049250313e-16;
626  IssmDouble meanvel = 3.170979198376458e-05; /*1000 m/yr*/
627 
628  /*Get Surface if required by one response*/
629  Input2* S_input = NULL;
630  for(int resp=0;resp<num_responses;resp++){
631  if(responses[resp]==SurfaceAverageVelMisfitEnum){
632  S_input = element->GetInput2(SurfaceAreaEnum); _assert_(S_input); break;
633  }
634  }
635 
636  /* Start looping on the number of gaussian points: */
637  Gauss* gauss=element->NewGaussTop(4);
638  for(int ig=gauss->begin();ig<gauss->end();ig++){
639  gauss->GaussPoint(ig);
640 
641  element->JacobianDeterminantTop(&Jdet,xyz_list_top,gauss);
642  element->NodalFunctions(basis, gauss);
643 
644  vx_input->GetInputValue(&vx,gauss);
645  vxobs_input->GetInputValue(&vxobs,gauss);
646  if(domaintype!=Domain2DverticalEnum){
647  vy_input->GetInputValue(&vy,gauss);
648  vyobs_input->GetInputValue(&vyobs,gauss);
649  }
650  /*Loop over all requested responses*/
651  for(int resp=0;resp<num_responses;resp++){
652  weights_input->GetInputValue(&weight,gauss,responses[resp]);
653 
654  switch(responses[resp]){
655  case SurfaceAbsVelMisfitEnum:
656  /*
657  * 1 [ 2 2 ]
658  * J = --- | (u - u ) + (v - v ) |
659  * 2 [ obs obs ]
660  *
661  * dJ
662  * DU = - -- = (u - u )
663  * du obs
664  */
665  for(i=0;i<numnodes;i++){
666  if(domaintype!=Domain2DverticalEnum){
667  dux=vxobs-vx;
668  duy=vyobs-vy;
669  pe->values[i*2+0]+=dux*weight*Jdet*gauss->weight*basis[i];
670  pe->values[i*2+1]+=duy*weight*Jdet*gauss->weight*basis[i];
671  }
672  else{
673  dux=vxobs-vx;
674  pe->values[i]+=dux*weight*Jdet*gauss->weight*basis[i];
675  }
676  }
677  break;
678  case SurfaceRelVelMisfitEnum:
679  /*
680  * 1 [ \bar{v}^2 2 \bar{v}^2 2 ]
681  * J = --- | ------------- (u - u ) + ------------- (v - v ) |
682  * 2 [ (u + eps)^2 obs (v + eps)^2 obs ]
683  * obs obs
684  *
685  * dJ \bar{v}^2
686  * DU = - -- = ------------- (u - u )
687  * du (u + eps)^2 obs
688  * obs
689  */
690  for(i=0;i<numnodes;i++){
691  if(domaintype!=Domain2DverticalEnum){
692  scalex=pow(meanvel/(vxobs+epsvel),2); if(vxobs==0)scalex=0;
693  scaley=pow(meanvel/(vyobs+epsvel),2); if(vyobs==0)scaley=0;
694  dux=scalex*(vxobs-vx);
695  duy=scaley*(vyobs-vy);
696  pe->values[i*2+0]+=dux*weight*Jdet*gauss->weight*basis[i];
697  pe->values[i*2+1]+=duy*weight*Jdet*gauss->weight*basis[i];
698  }
699  else{
700  scalex=pow(meanvel/(vxobs+epsvel),2); if(vxobs==0)scalex=0;
701  dux=scalex*(vxobs-vx);
702  pe->values[i]+=dux*weight*Jdet*gauss->weight*basis[i];
703  }
704  }
705  break;
706  case SurfaceLogVelMisfitEnum:
707  /*
708  * [ vel + eps ] 2
709  * J = 4 \bar{v}^2 | log ( ----------- ) |
710  * [ vel + eps ]
711  * obs
712  *
713  * dJ 2 * log(...)
714  * DU = - -- = - 4 \bar{v}^2 ------------- u
715  * du vel^2 + eps
716  *
717  */
718  for(i=0;i<numnodes;i++){
719  if(domaintype!=Domain2DverticalEnum){
720  velocity_mag =sqrt(pow(vx, 2)+pow(vy, 2))+epsvel;
721  obs_velocity_mag=sqrt(pow(vxobs,2)+pow(vyobs,2))+epsvel;
722  scale=-8*pow(meanvel,2)/pow(velocity_mag,2)*log(velocity_mag/obs_velocity_mag);
723  dux=scale*vx;
724  duy=scale*vy;
725  pe->values[i*2+0]+=dux*weight*Jdet*gauss->weight*basis[i];
726  pe->values[i*2+1]+=duy*weight*Jdet*gauss->weight*basis[i];
727  }
728  else{
729  velocity_mag =fabs(vx)+epsvel;
730  obs_velocity_mag=fabs(vxobs)+epsvel;
731  scale=-8*pow(meanvel,2)/pow(velocity_mag,2)*log(velocity_mag/obs_velocity_mag);
732  dux=scale*vx;
733  pe->values[i]+=dux*weight*Jdet*gauss->weight*basis[i];
734  }
735  }
736  break;
737  case SurfaceAverageVelMisfitEnum:
738  /*
739  * 1 2 2
740  * J = --- sqrt( (u - u ) + (v - v ) )
741  * S obs obs
742  *
743  * dJ 1 1
744  * DU = - -- = - --- ----------- * 2 (u - u )
745  * du S 2 sqrt(...) obs
746  */
747  S_input->GetInputValue(&S,gauss);
748  for(i=0;i<numnodes;i++){
749  if(domaintype!=Domain2DverticalEnum){
750  scale=1./(S*2*sqrt(pow(vx-vxobs,2)+pow(vy-vyobs,2))+epsvel);
751  dux=scale*(vxobs-vx);
752  duy=scale*(vyobs-vy);
753  pe->values[i*2+0]+=dux*weight*Jdet*gauss->weight*basis[i];
754  pe->values[i*2+1]+=duy*weight*Jdet*gauss->weight*basis[i];
755  }
756  else{
757  scale=1./(S*2*fabs(vx-vxobs)+epsvel);
758  dux=scale*(vxobs-vx);
759  pe->values[i]+=dux*weight*Jdet*gauss->weight*basis[i];
760  }
761  }
762  break;
763  case SurfaceLogVxVyMisfitEnum:
764  /*
765  * 1 [ |u| + eps 2 |v| + eps 2 ]
766  * J = --- \bar{v}^2 | log ( ----------- ) + log ( ----------- ) |
767  * 2 [ |u |+ eps |v |+ eps ]
768  * obs obs
769  * dJ 1 u 1
770  * DU = - -- = - \bar{v}^2 log(u...) --------- ---- ~ - \bar{v}^2 log(u...) ------
771  * du |u| + eps |u| u + eps
772  */
773  for(i=0;i<numnodes;i++){
774  if(domaintype!=Domain2DverticalEnum){
775  dux = - meanvel*meanvel * log((fabs(vx)+epsvel)/(fabs(vxobs)+epsvel)) / (vx+epsvel);
776  duy = - meanvel*meanvel * log((fabs(vy)+epsvel)/(fabs(vyobs)+epsvel)) / (vy+epsvel);
777  pe->values[i*2+0]+=dux*weight*Jdet*gauss->weight*basis[i];
778  pe->values[i*2+1]+=duy*weight*Jdet*gauss->weight*basis[i];
779  }
780  else{
781  dux = - meanvel*meanvel * log((fabs(vx)+epsvel)/(fabs(vxobs)+epsvel)) / (vx+epsvel);
782  pe->values[i]+=dux*weight*Jdet*gauss->weight*basis[i];
783  }
784  }
785  break;
786  case DragCoefficientAbsGradientEnum:
787  /*Nothing in P vector*/
788  break;
789  case ThicknessAbsGradientEnum:
790  /*Nothing in P vector*/
791  break;
792  case ThicknessAlongGradientEnum:
793  /*Nothing in P vector*/
794  break;
795  case ThicknessAcrossGradientEnum:
796  /*Nothing in P vector*/
797  break;
798  case RheologyBbarAbsGradientEnum:
799  /*Nothing in P vector*/
800  break;
801  case RheologyBAbsGradientEnum:
802  /*Nothing in P vector*/
803  break;
804  case RheologyBInitialguessMisfitEnum:
805  /*Nothing in P vector*/
806  break;
807  default:
808  _error_("response " << EnumToStringx(responses[resp]) << " not supported yet");
809  }
810  }
811  }
812 
813  /*Transform coordinate system*/
814  if(domaintype!=Domain2DverticalEnum) element->TransformLoadVectorCoord(pe,XYEnum);
815 
816  /*Clean up and return*/
817  xDelete<int>(responses);
818  xDelete<IssmDouble>(xyz_list_top);
819  xDelete<IssmDouble>(basis);
820  delete gauss;
821  return pe;
822 
823 }/*}}}*/

◆ CreatePVectorSSA()

ElementVector * AdjointHorizAnalysis::CreatePVectorSSA ( Element *  element)

Definition at line 824 of file AdjointHorizAnalysis.cpp.

824  {/*{{{*/
825 
826  /*Intermediaries*/
827  int domaintype;
828  Element* basalelement;
829 
830  /* Check if ice in element */
831  if(!element->IsIceInElement()) return NULL;
832 
833  /*Get basal element*/
834  element->FindParam(&domaintype,DomainTypeEnum);
835  switch(domaintype){
836  case Domain2DhorizontalEnum:
837  basalelement = element;
838  break;
839  case Domain2DverticalEnum:
840  if(!element->IsOnBase()) return NULL;
841  basalelement = element->SpawnBasalElement();
842  break;
843  case Domain3DEnum:
844  if(!element->IsOnBase()) return NULL;
845  basalelement = element->SpawnBasalElement();
846  break;
847  default: _error_("mesh "<<EnumToStringx(domaintype)<<" not supported yet");
848  }
849 
850  /*Intermediaries */
851  int num_responses,i;
852  IssmDouble Jdet,obs_velocity_mag,velocity_mag;
853  IssmDouble vx,vy,vxobs,vyobs,dux,duy,weight;
854  IssmDouble scalex,scaley,scale,S;
855  int *responses = NULL;
856  IssmDouble *xyz_list = NULL;
857 
858  /*Fetch number of nodes and dof for this finite element*/
859  int numnodes = basalelement->GetNumberOfNodes();
860 
861  /*Initialize Element vector and vectors*/
862  ElementVector* pe = basalelement->NewElementVector(SSAApproximationEnum);
863  IssmDouble* basis = xNew<IssmDouble>(numnodes);
864 
865  /*Retrieve all inputs and parameters*/
866  basalelement->GetVerticesCoordinates(&xyz_list);
867  basalelement->FindParam(&num_responses,InversionNumCostFunctionsEnum);
868  basalelement->FindParam(&responses,NULL,InversionCostFunctionsEnum);
869  DatasetInput2* weights_input = basalelement->GetDatasetInput2(InversionCostFunctionsCoefficientsEnum); _assert_(weights_input);
870  Input2* vx_input = basalelement->GetInput2(VxEnum); _assert_(vx_input);
871  Input2* vxobs_input = basalelement->GetInput2(InversionVxObsEnum); _assert_(vxobs_input);
872  Input2* vy_input=NULL;
873  Input2* vyobs_input=NULL;
874  if(domaintype!=Domain2DverticalEnum){
875  vy_input = basalelement->GetInput2(VyEnum); _assert_(vy_input);
876  vyobs_input = basalelement->GetInput2(InversionVyObsEnum); _assert_(vyobs_input);
877  }
878  IssmDouble epsvel = 2.220446049250313e-16;
879  IssmDouble meanvel = 3.170979198376458e-05; /*1000 m/yr*/
880 
881  /*Get Surface if required by one response*/
882  Input2* S_input = NULL;
883  for(int resp=0;resp<num_responses;resp++){
884  if(responses[resp]==SurfaceAverageVelMisfitEnum){
885  S_input = element->GetInput2(SurfaceAreaEnum); _assert_(S_input); break;
886  }
887  }
888 
889  /* Start looping on the number of gaussian points: */
890  Gauss* gauss=basalelement->NewGauss(4);
891  for(int ig=gauss->begin();ig<gauss->end();ig++){
892  gauss->GaussPoint(ig);
893 
894  basalelement->JacobianDeterminant(&Jdet,xyz_list,gauss);
895  basalelement->NodalFunctions(basis, gauss);
896 
897  vx_input->GetInputValue(&vx,gauss);
898  vxobs_input->GetInputValue(&vxobs,gauss);
899  if(domaintype!=Domain2DverticalEnum){
900  vy_input->GetInputValue(&vy,gauss);
901  vyobs_input->GetInputValue(&vyobs,gauss);
902  }
903  /*Loop over all requested responses*/
904  for(int resp=0;resp<num_responses;resp++){
905  weights_input->GetInputValue(&weight,gauss,responses[resp]);
906 
907  switch(responses[resp]){
908  case SurfaceAbsVelMisfitEnum:
909  /*
910  * 1 [ 2 2 ]
911  * J = --- | (u - u ) + (v - v ) |
912  * 2 [ obs obs ]
913  *
914  * dJ
915  * DU = - -- = (u - u )
916  * du obs
917  */
918  for(i=0;i<numnodes;i++){
919  if(domaintype!=Domain2DverticalEnum){
920  dux=vxobs-vx;
921  duy=vyobs-vy;
922  pe->values[i*2+0]+=dux*weight*Jdet*gauss->weight*basis[i];
923  pe->values[i*2+1]+=duy*weight*Jdet*gauss->weight*basis[i];
924  }
925  else {
926  dux=vxobs-vx;
927  pe->values[i]+=dux*weight*Jdet*gauss->weight*basis[i];
928  }
929  }
930  break;
931  case SurfaceRelVelMisfitEnum:
932  /*
933  * 1 [ \bar{v}^2 2 \bar{v}^2 2 ]
934  * J = --- | ------------- (u - u ) + ------------- (v - v ) |
935  * 2 [ (u + eps)^2 obs (v + eps)^2 obs ]
936  * obs obs
937  *
938  * dJ \bar{v}^2
939  * DU = - -- = ------------- (u - u )
940  * du (u + eps)^2 obs
941  * obs
942  */
943  for(i=0;i<numnodes;i++){
944  if(domaintype!=Domain2DverticalEnum){
945  scalex=pow(meanvel/(vxobs+epsvel),2); if(vxobs==0)scalex=0;
946  scaley=pow(meanvel/(vyobs+epsvel),2); if(vyobs==0)scaley=0;
947  dux=scalex*(vxobs-vx);
948  duy=scaley*(vyobs-vy);
949  pe->values[i*2+0]+=dux*weight*Jdet*gauss->weight*basis[i];
950  pe->values[i*2+1]+=duy*weight*Jdet*gauss->weight*basis[i];
951  }
952  else{
953  scalex=pow(meanvel/(vxobs+epsvel),2); if(vxobs==0)scalex=0;
954  dux=scalex*(vxobs-vx);
955  pe->values[i]+=dux*weight*Jdet*gauss->weight*basis[i];
956  }
957  }
958  break;
959  case SurfaceLogVelMisfitEnum:
960  /*
961  * [ vel + eps ] 2
962  * J = 4 \bar{v}^2 | log ( ----------- ) |
963  * [ vel + eps ]
964  * obs
965  *
966  * dJ 2 * log(...)
967  * DU = - -- = - 4 \bar{v}^2 ------------- u
968  * du vel^2 + eps
969  *
970  */
971  for(i=0;i<numnodes;i++){
972  if(domaintype!=Domain2DverticalEnum){
973  velocity_mag =sqrt(pow(vx, 2)+pow(vy, 2))+epsvel;
974  obs_velocity_mag=sqrt(pow(vxobs,2)+pow(vyobs,2))+epsvel;
975  scale=-8*pow(meanvel,2)/pow(velocity_mag,2)*log(velocity_mag/obs_velocity_mag);
976  dux=scale*vx;
977  duy=scale*vy;
978  pe->values[i*2+0]+=dux*weight*Jdet*gauss->weight*basis[i];
979  pe->values[i*2+1]+=duy*weight*Jdet*gauss->weight*basis[i];
980  }
981  else{
982  velocity_mag =fabs(vx)+epsvel;
983  obs_velocity_mag=fabs(vxobs)+epsvel;
984  scale=-8*pow(meanvel,2)/pow(velocity_mag,2)*log(velocity_mag/obs_velocity_mag);
985  dux=scale*vx;
986  pe->values[i]+=dux*weight*Jdet*gauss->weight*basis[i];
987  }
988  }
989  break;
990  case SurfaceAverageVelMisfitEnum:
991  /*
992  * 1 2 2
993  * J = --- sqrt( (u - u ) + (v - v ) )
994  * S obs obs
995  *
996  * dJ 1 1
997  * DU = - -- = - --- ----------- * 2 (u - u )
998  * du S 2 sqrt(...) obs
999  */
1000  S_input->GetInputValue(&S,gauss);
1001  for(i=0;i<numnodes;i++){
1002  if(domaintype!=Domain2DverticalEnum){
1003  scale=1./(S*2*sqrt(pow(vx-vxobs,2)+pow(vy-vyobs,2))+epsvel);
1004  dux=scale*(vxobs-vx);
1005  duy=scale*(vyobs-vy);
1006  pe->values[i*2+0]+=dux*weight*Jdet*gauss->weight*basis[i];
1007  pe->values[i*2+1]+=duy*weight*Jdet*gauss->weight*basis[i];
1008  }
1009  else{
1010  scale=1./(S*2*fabs(vx-vxobs)+epsvel);
1011  dux=scale*(vxobs-vx);
1012  pe->values[i]+=dux*weight*Jdet*gauss->weight*basis[i];
1013  }
1014  }
1015  break;
1016  case SurfaceLogVxVyMisfitEnum:
1017  /*
1018  * 1 [ |u| + eps 2 |v| + eps 2 ]
1019  * J = --- \bar{v}^2 | log ( ----------- ) + log ( ----------- ) |
1020  * 2 [ |u |+ eps |v |+ eps ]
1021  * obs obs
1022  * dJ 1 u 1
1023  * DU = - -- = - \bar{v}^2 log(u...) --------- ---- ~ - \bar{v}^2 log(u...) ------
1024  * du |u| + eps |u| u + eps
1025  */
1026  for(i=0;i<numnodes;i++){
1027  if(domaintype!=Domain2DverticalEnum){
1028  dux = - meanvel*meanvel * log((fabs(vx)+epsvel)/(fabs(vxobs)+epsvel)) / (vx+epsvel);
1029  duy = - meanvel*meanvel * log((fabs(vy)+epsvel)/(fabs(vyobs)+epsvel)) / (vy+epsvel);
1030  pe->values[i*2+0]+=dux*weight*Jdet*gauss->weight*basis[i];
1031  pe->values[i*2+1]+=duy*weight*Jdet*gauss->weight*basis[i];
1032  }
1033  else{
1034  dux = - meanvel*meanvel * log((fabs(vx)+epsvel)/(fabs(vxobs)+epsvel)) / (vx+epsvel);
1035  pe->values[i]+=dux*weight*Jdet*gauss->weight*basis[i];
1036  }
1037  }
1038  break;
1039  case DragCoefficientAbsGradientEnum:
1040  /*Nothing in P vector*/
1041  break;
1042  case ThicknessAbsGradientEnum:
1043  /*Nothing in P vector*/
1044  break;
1045  case ThicknessAlongGradientEnum:
1046  /*Nothing in P vector*/
1047  break;
1048  case ThicknessAcrossGradientEnum:
1049  /*Nothing in P vector*/
1050  break;
1051  case RheologyBbarAbsGradientEnum:
1052  /*Nothing in P vector*/
1053  break;
1054  case RheologyBAbsGradientEnum:
1055  /*Nothing in P vector*/
1056  break;
1057  case RheologyBInitialguessMisfitEnum:
1058  /*Nothing in P vector*/
1059  break;
1060  default:
1061  _error_("response " << EnumToStringx(responses[resp]) << " not supported yet");
1062  }
1063  }
1064  }
1065 
1066  /*Transform coordinate system*/
1067  if(domaintype!=Domain2DverticalEnum) basalelement->TransformLoadVectorCoord(pe,XYEnum);
1068 
1069  /*Clean up and return*/
1070  xDelete<int>(responses);
1071  xDelete<IssmDouble>(xyz_list);
1072  xDelete<IssmDouble>(basis);
1073  if(domaintype!=Domain2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;};
1074  delete gauss;
1075  return pe;
1076 }/*}}}*/

◆ GetSolutionFromInputs()

void AdjointHorizAnalysis::GetSolutionFromInputs ( Vector< IssmDouble > *  solution,
Element *  element 
)
virtual

Implements Analysis.

Definition at line 1077 of file AdjointHorizAnalysis.cpp.

1077  {/*{{{*/
1078  _error_("not implemented yet");
1079 }/*}}}*/

◆ GradientJ()

void AdjointHorizAnalysis::GradientJ ( Vector< IssmDouble > *  gradient,
Element *  element,
int  control_type,
int  control_index 
)
virtual

Implements Analysis.

Definition at line 1080 of file AdjointHorizAnalysis.cpp.

1080  {/*{{{*/
1081  /*The gradient of the cost function is calculated in 2 parts.
1082  *
1083  * dJ \partial J \partial lambda^T(KU-F)
1084  * -- = ---------- + ------------------------
1085  * dk \partial k \parial k
1086  *
1087  * */
1088 
1089  /*If on water, grad = 0: */
1090  if(!element->IsIceInElement()) return;
1091 
1092  /*Get Approximation*/
1093  int approximation;
1094  element->GetInputValue(&approximation,ApproximationEnum);
1095 
1096  /*Get list of cost functions*/
1097  int *responses = NULL;
1098  int num_responses,resp;
1099  element->FindParam(&num_responses,InversionNumCostFunctionsEnum);
1100  element->FindParam(&responses,NULL,InversionCostFunctionsEnum);
1101 
1102  /*Check that control_type is supported*/
1103  if(control_type!=MaterialsRheologyBbarEnum &&
1104  control_type!=FrictionCoefficientEnum &&
1105  control_type!=FrictionAsEnum &&
1106  control_type!=DamageDbarEnum &&
1107  control_type!=MaterialsRheologyBEnum){
1108  _error_("Control "<<EnumToStringx(control_type)<<" not supported");
1109  }
1110 
1111  /*Deal with first part (partial derivative a J with respect to k)*/
1112  for(resp=0;resp<num_responses;resp++) switch(responses[resp]){
1113  case SurfaceAbsVelMisfitEnum: /*Nothing, \partial J/\partial k = 0*/ break;
1114  case SurfaceRelVelMisfitEnum: /*Nothing, \partial J/\partial k = 0*/ break;
1115  case SurfaceLogVelMisfitEnum: /*Nothing, \partial J/\partial k = 0*/ break;
1116  case SurfaceLogVxVyMisfitEnum: /*Nothing, \partial J/\partial k = 0*/ break;
1117  case SurfaceAverageVelMisfitEnum: /*Nothing, \partial J/\partial k = 0*/ break;
1118  case DragCoefficientAbsGradientEnum: GradientJDragGradient(element,gradient,control_index); break;
1119  case RheologyBbarAbsGradientEnum: GradientJBbarGradient(element,gradient,control_index); break;
1120  case RheologyBAbsGradientEnum: GradientJBGradient(element,gradient,control_index); break;
1121  case RheologyBInitialguessMisfitEnum: GradientJBinitial(element,gradient,control_index); break;
1122  default: _error_("response " << EnumToStringx(responses[resp]) << " not supported yet");
1123  }
1124 
1125  /*Deal with second term*/
1126  switch(control_type){
1127  case FrictionCoefficientEnum:
1128  switch(approximation){
1129  case SSAApproximationEnum: GradientJDragSSA(element,gradient,control_index); break;
1130  case L1L2ApproximationEnum:GradientJDragL1L2(element,gradient,control_index); break;
1131  case HOApproximationEnum: GradientJDragHO( element,gradient,control_index); break;
1132  case FSApproximationEnum: GradientJDragFS( element,gradient,control_index); break;
1133  case NoneApproximationEnum: /*Gradient is 0*/ break;
1134  default: _error_("approximation " << EnumToStringx(approximation) << " not supported yet");
1135  }
1136  break;
1137  case FrictionAsEnum:
1138  switch(approximation){
1139  case SSAApproximationEnum: GradientJDragHydroSSA(element,gradient,control_index); break;
1140  case L1L2ApproximationEnum:GradientJDragHydroL1L2(element,gradient,control_index); break;
1141  case HOApproximationEnum: GradientJDragHydroHO( element,gradient,control_index); break;
1142  case FSApproximationEnum: GradientJDragHydroFS( element,gradient,control_index); break;
1143  case NoneApproximationEnum: /*Gradient is 0*/ break;
1144  default: _error_("approximation " << EnumToStringx(approximation) << " not supported yet");
1145  }
1146  break;
1147  case MaterialsRheologyBbarEnum:
1148  switch(approximation){
1149  case SSAApproximationEnum: GradientJBbarSSA(element,gradient,control_index); break;
1150  case L1L2ApproximationEnum:GradientJBbarL1L2(element,gradient,control_index); break;
1151  case HOApproximationEnum: GradientJBbarHO( element,gradient,control_index); break;
1152  case FSApproximationEnum: GradientJBbarFS( element,gradient,control_index); break;
1153  case NoneApproximationEnum: /*Gradient is 0*/ break;
1154  default: _error_("approximation " << EnumToStringx(approximation) << " not supported yet");
1155  }
1156  break;
1157  case MaterialsRheologyBEnum:
1158  switch(approximation){
1159  case SSAApproximationEnum: GradientJBSSA(element,gradient,control_index); break;
1160  case HOApproximationEnum: GradientJBHO( element,gradient,control_index); break;
1161  case FSApproximationEnum: GradientJBFS( element,gradient,control_index); break;
1162  case NoneApproximationEnum: /*Gradient is 0*/ break;
1163  default: _error_("approximation " << EnumToStringx(approximation) << " not supported yet");
1164  }
1165  break;
1166  case DamageDbarEnum:
1167  switch(approximation){
1168  case SSAApproximationEnum: GradientJDSSA(element,gradient,control_index); break;
1169  case NoneApproximationEnum: /*Gradient is 0*/ break;
1170  default: _error_("approximation " << EnumToStringx(approximation) << " not supported yet");
1171  }
1172  break;
1173  default: _error_("control type not supported yet: " << EnumToStringx(control_type));
1174  }
1175 
1176  /*Clean up and return*/
1177  xDelete<int>(responses);
1178 
1179 }/*}}}*/

◆ GradientJBbarFS()

void AdjointHorizAnalysis::GradientJBbarFS ( Element *  element,
Vector< IssmDouble > *  gradient,
int  control_index 
)

Definition at line 1180 of file AdjointHorizAnalysis.cpp.

1180  {/*{{{*/
1181  /*WARNING: We use SSA as an estimate for now*/
1182  this->GradientJBbarSSA(element,gradient,control_index);
1183 }/*}}}*/

◆ GradientJBbarGradient()

void AdjointHorizAnalysis::GradientJBbarGradient ( Element *  element,
Vector< IssmDouble > *  gradient,
int  control_index 
)

Definition at line 1184 of file AdjointHorizAnalysis.cpp.

1184  {/*{{{*/
1185 
1186  /*Intermediaries*/
1187  int domaintype;
1188  Element* basalelement;
1189 
1190  /*Get basal element*/
1191  element->FindParam(&domaintype,DomainTypeEnum);
1192  switch(domaintype){
1193  case Domain2DhorizontalEnum:
1194  basalelement = element;
1195  break;
1196  case Domain2DverticalEnum:
1197  if(!element->IsOnBase()) return;
1198  basalelement = element->SpawnBasalElement();
1199  break;
1200  case Domain3DEnum:
1201  if(!element->IsOnBase()) return;
1202  basalelement = element->SpawnBasalElement();
1203  break;
1204  default: _error_("mesh "<<EnumToStringx(domaintype)<<" not supported yet");
1205  }
1206 
1207  /*Intermediaries*/
1208  IssmDouble Jdet,weight;
1209  IssmDouble dk[3];
1210  IssmDouble *xyz_list= NULL;
1211 
1212  /*Fetch number of vertices for this finite element*/
1213  int numvertices = basalelement->GetNumberOfVertices();
1214 
1215  /*Initialize some vectors*/
1216  IssmDouble* dbasis = xNew<IssmDouble>(2*numvertices);
1217  IssmDouble* ge = xNewZeroInit<IssmDouble>(numvertices);
1218  int* vertexpidlist = xNew<int>(numvertices);
1219 
1220  /*Retrieve all inputs we will be needing: */
1221  basalelement->GetVerticesCoordinates(&xyz_list);
1222  basalelement->GradientIndexing(&vertexpidlist[0],control_index);
1223  Input2* rheologyb_input = basalelement->GetInput2(MaterialsRheologyBbarEnum); _assert_(rheologyb_input);
1224  DatasetInput2* weights_input = basalelement->GetDatasetInput2(InversionCostFunctionsCoefficientsEnum); _assert_(weights_input);
1225 
1226  /* Start looping on the number of gaussian points: */
1227  Gauss* gauss=basalelement->NewGauss(2);
1228  for(int ig=gauss->begin();ig<gauss->end();ig++){
1229  gauss->GaussPoint(ig);
1230 
1231  basalelement->JacobianDeterminant(&Jdet,xyz_list,gauss);
1232  basalelement->NodalFunctionsP1Derivatives(dbasis,xyz_list,gauss);
1233  weights_input->GetInputValue(&weight,gauss,RheologyBbarAbsGradientEnum);
1234 
1235  /*Build alpha_complement_list: */
1236  rheologyb_input->GetInputDerivativeValue(&dk[0],xyz_list,gauss);
1237 
1238  /*Build gradje_g_gaussian vector (actually -dJ/ddrag): */
1239  for(int i=0;i<numvertices;i++){
1240  if(domaintype!=Domain2DverticalEnum){
1241  ge[i]+=-weight*Jdet*gauss->weight*(dbasis[0*numvertices+i]*dk[0]+dbasis[1*numvertices+i]*dk[1]);
1242  }
1243  else{
1244  ge[i]+=-weight*Jdet*gauss->weight*dbasis[0*numvertices+i]*dk[0];
1245  }
1246  _assert_(!xIsNan<IssmDouble>(ge[i]));
1247  }
1248  }
1249  gradient->SetValues(numvertices,vertexpidlist,ge,ADD_VAL);
1250 
1251  /*Clean up and return*/
1252  xDelete<IssmDouble>(xyz_list);
1253  xDelete<IssmDouble>(dbasis);
1254  xDelete<IssmDouble>(ge);
1255  xDelete<int>(vertexpidlist);
1256  delete gauss;
1257  if(domaintype!=Domain2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;};
1258 }/*}}}*/

◆ GradientJBinitial()

void AdjointHorizAnalysis::GradientJBinitial ( Element *  element,
Vector< IssmDouble > *  gradient,
int  control_index 
)

Definition at line 1423 of file AdjointHorizAnalysis.cpp.

1423  {/*{{{*/
1424 
1425  /*Intermediaries*/
1426  int domaintype;
1427 
1428  /*Get basal element*/
1429  element->FindParam(&domaintype,DomainTypeEnum);
1430  switch(domaintype){
1431  case Domain2DhorizontalEnum:
1432  break;
1433  case Domain2DverticalEnum:
1434  break;
1435  case Domain3DEnum:
1436  break;
1437  default: _error_("mesh "<<EnumToStringx(domaintype)<<" not supported yet");
1438  }
1439 
1440  /*Intermediaries*/
1441  IssmDouble Jdet,weight;
1442  IssmDouble B,B0;
1443  IssmDouble *xyz_list= NULL;
1444 
1445  /*Fetch number of vertices for this finite element*/
1446  int numvertices = element->GetNumberOfVertices();
1447 
1448  /*Initialize some vectors*/
1449  IssmDouble* basis = xNew<IssmDouble>(numvertices);
1450  IssmDouble* ge = xNewZeroInit<IssmDouble>(numvertices);
1451  int* vertexpidlist = xNew<int>(numvertices);
1452 
1453  /*Retrieve all inputs we will be needing: */
1454  element->GetVerticesCoordinates(&xyz_list);
1455  element->GradientIndexing(&vertexpidlist[0],control_index);
1456  Input2* rheology_input = element->GetInput2(MaterialsRheologyBbarEnum); _assert_(rheology_input);
1457  Input2* rheology0_input = element->GetInput2(RheologyBInitialguessEnum); _assert_(rheology0_input);
1458  DatasetInput2* weights_input = element->GetDatasetInput2(InversionCostFunctionsCoefficientsEnum); _assert_(weights_input);
1459 
1460  /* Start looping on the number of gaussian points: */
1461  Gauss* gauss=element->NewGauss(2);
1462  for(int ig=gauss->begin();ig<gauss->end();ig++){
1463  gauss->GaussPoint(ig);
1464  element->JacobianDeterminant(&Jdet,xyz_list,gauss);
1465  element->NodalFunctionsP1(basis,gauss);
1466  weights_input->GetInputValue(&weight,gauss,RheologyBInitialguessMisfitEnum);
1467 
1468  /*Build alpha_complement_list: */
1469  rheology_input->GetInputValue(&B,gauss);
1470  rheology0_input->GetInputValue(&B0,gauss);
1471 
1472  /*Build gradje_g_gaussian vector (actually -dJ/ddrag): */
1473  for(int i=0;i<numvertices;i++){
1474  ge[i]+=-weight*Jdet*gauss->weight*basis[i]*(B-B0);
1475  _assert_(!xIsNan<IssmDouble>(ge[i]));
1476  }
1477  }
1478  gradient->SetValues(numvertices,vertexpidlist,ge,ADD_VAL);
1479 
1480  /*Clean up and return*/
1481  xDelete<IssmDouble>(xyz_list);
1482  xDelete<IssmDouble>(basis);
1483  xDelete<IssmDouble>(ge);
1484  xDelete<int>(vertexpidlist);
1485  delete gauss;
1486 }/*}}}*/

◆ GradientJBbarL1L2()

void AdjointHorizAnalysis::GradientJBbarL1L2 ( Element *  element,
Vector< IssmDouble > *  gradient,
int  control_index 
)

Definition at line 1259 of file AdjointHorizAnalysis.cpp.

1259  {/*{{{*/
1260 
1261  /*Same as SSA*/
1262  return this->GradientJBbarSSA(element,gradient,control_index);
1263 }/*}}}*/

◆ GradientJBbarHO()

void AdjointHorizAnalysis::GradientJBbarHO ( Element *  element,
Vector< IssmDouble > *  gradient,
int  control_index 
)

Definition at line 1264 of file AdjointHorizAnalysis.cpp.

1264  {/*{{{*/
1265 
1266  /*WARNING: We use SSA as an estimate for now*/
1267  this->GradientJBbarSSA(element,gradient,control_index);
1268 }/*}}}*/

◆ GradientJBbarSSA()

void AdjointHorizAnalysis::GradientJBbarSSA ( Element *  element,
Vector< IssmDouble > *  gradient,
int  control_index 
)

Definition at line 1269 of file AdjointHorizAnalysis.cpp.

1269  {/*{{{*/
1270 
1271  /*Intermediaries*/
1272  int domaintype,dim;
1273  Element* basalelement;
1274 
1275  /*Get basal element*/
1276  element->FindParam(&domaintype,DomainTypeEnum);
1277  switch(domaintype){
1278  case Domain2DhorizontalEnum:
1279  basalelement = element;
1280  dim = 2;
1281  break;
1282  case Domain2DverticalEnum:
1283  if(!element->IsOnBase()) return;
1284  basalelement = element->SpawnBasalElement();
1285  dim = 1;
1286  break;
1287  case Domain3DEnum:
1288  if(!element->IsOnBase()) return;
1289  basalelement = element->SpawnBasalElement();
1290  dim = 2;
1291  break;
1292  default: _error_("mesh "<<EnumToStringx(domaintype)<<" not supported yet");
1293  }
1294 
1295  /*Intermediaries*/
1296  IssmDouble Jdet,weight;
1297  IssmDouble thickness,dmudB;
1298  IssmDouble dvx[3],dvy[3],dadjx[3],dadjy[3];
1299  IssmDouble *xyz_list= NULL;
1300 
1301  /*Fetch number of vertices for this finite element*/
1302  int numvertices = basalelement->GetNumberOfVertices();
1303 
1304  /*Initialize some vectors*/
1305  IssmDouble* basis = xNew<IssmDouble>(numvertices);
1306  IssmDouble* ge = xNewZeroInit<IssmDouble>(numvertices);
1307  int* vertexpidlist = xNew<int>(numvertices);
1308 
1309  /*Retrieve all inputs we will be needing: */
1310  basalelement->GetVerticesCoordinates(&xyz_list);
1311  basalelement->GradientIndexing(&vertexpidlist[0],control_index);
1312  Input2* thickness_input = basalelement->GetInput2(ThicknessEnum); _assert_(thickness_input);
1313  Input2* vx_input = basalelement->GetInput2(VxEnum); _assert_(vx_input);
1314  Input2* vy_input = basalelement->GetInput2(VyEnum); _assert_(vy_input);
1315  Input2* adjointx_input = basalelement->GetInput2(AdjointxEnum); _assert_(adjointx_input);
1316  Input2* adjointy_input = basalelement->GetInput2(AdjointyEnum); _assert_(adjointy_input);
1317  Input2* rheologyb_input = basalelement->GetInput2(MaterialsRheologyBbarEnum); _assert_(rheologyb_input);
1318 
1319  /* Start looping on the number of gaussian points: */
1320  Gauss* gauss=basalelement->NewGauss(4);
1321  for(int ig=gauss->begin();ig<gauss->end();ig++){
1322  gauss->GaussPoint(ig);
1323 
1324  thickness_input->GetInputValue(&thickness,gauss);
1325  vx_input->GetInputDerivativeValue(&dvx[0],xyz_list,gauss);
1326  vy_input->GetInputDerivativeValue(&dvy[0],xyz_list,gauss);
1327  adjointx_input->GetInputDerivativeValue(&dadjx[0],xyz_list,gauss);
1328  adjointy_input->GetInputDerivativeValue(&dadjy[0],xyz_list,gauss);
1329 
1330  basalelement->dViscositydBSSA(&dmudB,dim,xyz_list,gauss,vx_input,vy_input);
1331 
1332  basalelement->JacobianDeterminant(&Jdet,xyz_list,gauss);
1333  basalelement->NodalFunctionsP1(basis,gauss);
1334 
1335  /*Build gradient vector (actually -dJ/dB): */
1336  for(int i=0;i<numvertices;i++){
1337  ge[i]+=-dmudB*thickness*(
1338  (2*dvx[0]+dvy[1])*2*dadjx[0]+(dvx[1]+dvy[0])*(dadjx[1]+dadjy[0])+(2*dvy[1]+dvx[0])*2*dadjy[1]
1339  )*Jdet*gauss->weight*basis[i];
1340  _assert_(!xIsNan<IssmDouble>(ge[i]));
1341  }
1342  }
1343  gradient->SetValues(numvertices,vertexpidlist,ge,ADD_VAL);
1344 
1345  /*Clean up and return*/
1346  xDelete<IssmDouble>(xyz_list);
1347  xDelete<IssmDouble>(basis);
1348  xDelete<IssmDouble>(ge);
1349  xDelete<int>(vertexpidlist);
1350  delete gauss;
1351  if(domaintype!=Domain2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;};
1352 }/*}}}*/

◆ GradientJBFS()

void AdjointHorizAnalysis::GradientJBFS ( Element *  element,
Vector< IssmDouble > *  gradient,
int  control_index 
)

Definition at line 1353 of file AdjointHorizAnalysis.cpp.

1353  {/*{{{*/
1354  /*WARNING: We use HO as an estimate for now*/
1355  this->GradientJBHO(element,gradient,control_index);
1356 }/*}}}*/

◆ GradientJBGradient()

void AdjointHorizAnalysis::GradientJBGradient ( Element *  element,
Vector< IssmDouble > *  gradient,
int  control_index 
)

Definition at line 1357 of file AdjointHorizAnalysis.cpp.

1357  {/*{{{*/
1358 
1359  /*Intermediaries*/
1360  int domaintype;
1361 
1362  /*Get basal element*/
1363  element->FindParam(&domaintype,DomainTypeEnum);
1364  switch(domaintype){
1365  case Domain2DhorizontalEnum:
1366  break;
1367  case Domain2DverticalEnum:
1368  break;
1369  case Domain3DEnum:
1370  break;
1371  default: _error_("mesh "<<EnumToStringx(domaintype)<<" not supported yet");
1372  }
1373 
1374  /*Intermediaries*/
1375  IssmDouble Jdet,weight;
1376  IssmDouble dk[3];
1377  IssmDouble *xyz_list= NULL;
1378 
1379  /*Fetch number of vertices for this finite element*/
1380  int numvertices = element->GetNumberOfVertices();
1381 
1382  /*Initialize some vectors*/
1383  IssmDouble* dbasis = xNew<IssmDouble>(3*numvertices);
1384  IssmDouble* ge = xNewZeroInit<IssmDouble>(numvertices);
1385  int* vertexpidlist = xNew<int>(numvertices);
1386 
1387  /*Retrieve all inputs we will be needing: */
1388  element->GetVerticesCoordinates(&xyz_list);
1389  element->GradientIndexing(&vertexpidlist[0],control_index);
1390  Input2* rheology_input = element->GetInput2(MaterialsRheologyBEnum); _assert_(rheology_input);
1391  DatasetInput2* weights_input = element->GetDatasetInput2(InversionCostFunctionsCoefficientsEnum); _assert_(weights_input);
1392  /* Start looping on the number of gaussian points: */
1393  Gauss* gauss=element->NewGauss(2);
1394  for(int ig=gauss->begin();ig<gauss->end();ig++){
1395  gauss->GaussPoint(ig);
1396  element->JacobianDeterminant(&Jdet,xyz_list,gauss);
1397  element->NodalFunctionsP1Derivatives(dbasis,xyz_list,gauss);
1398  weights_input->GetInputValue(&weight,gauss,RheologyBAbsGradientEnum);
1399 
1400  /*Build alpha_complement_list: */
1401  rheology_input->GetInputDerivativeValue(&dk[0],xyz_list,gauss);
1402 
1403  /*Build gradje_g_gaussian vector (actually -dJ/ddrag): */
1404  for(int i=0;i<numvertices;i++){
1405  if(domaintype!=Domain2DverticalEnum){
1406  ge[i]+=-weight*Jdet*gauss->weight*(dbasis[0*numvertices+i]*dk[0]+dbasis[1*numvertices+i]*dk[1]);
1407  }
1408  else{
1409  ge[i]+=-weight*Jdet*gauss->weight*(dbasis[0*numvertices+i]*dk[0]+dbasis[1*numvertices+i]*dk[1]+dbasis[2*numvertices+i]*dk[2]);
1410  }
1411  _assert_(!xIsNan<IssmDouble>(ge[i]));
1412  }
1413  }
1414  gradient->SetValues(numvertices,vertexpidlist,ge,ADD_VAL);
1415 
1416  /*Clean up and return*/
1417  xDelete<IssmDouble>(xyz_list);
1418  xDelete<IssmDouble>(dbasis);
1419  xDelete<IssmDouble>(ge);
1420  xDelete<int>(vertexpidlist);
1421  delete gauss;
1422 }/*}}}*/

◆ GradientJBHO()

void AdjointHorizAnalysis::GradientJBHO ( Element *  element,
Vector< IssmDouble > *  gradient,
int  control_index 
)

Definition at line 1487 of file AdjointHorizAnalysis.cpp.

1487  {/*{{{*/
1488  /*Intermediaries*/
1489  int domaintype,dim;
1490 
1491  /*Get domaintype*/
1492  element->FindParam(&domaintype,DomainTypeEnum);
1493 
1494  /*Intermediaries*/
1495  IssmDouble Jdet,weight;
1496  IssmDouble thickness,dmudB;
1497  IssmDouble dvx[3],dvy[3],dadjx[3],dadjy[3];
1498  IssmDouble *xyz_list= NULL;
1499 
1500  /*Fetch number of vertices for this finite element*/
1501  int numvertices = element->GetNumberOfVertices();
1502 
1503  /*Initialize some vectors*/
1504  IssmDouble* basis = xNew<IssmDouble>(numvertices);
1505  IssmDouble* ge = xNewZeroInit<IssmDouble>(numvertices);
1506  int* vertexpidlist = xNew<int>(numvertices);
1507 
1508  /*Retrieve all inputs we will be needing: */
1509  element->GetVerticesCoordinates(&xyz_list);
1510  element->GradientIndexing(&vertexpidlist[0],control_index);
1511  Input2* thickness_input = element->GetInput2(ThicknessEnum); _assert_(thickness_input);
1512  Input2* vx_input = element->GetInput2(VxEnum); _assert_(vx_input);
1513  Input2* vy_input = NULL;
1514  Input2* adjointx_input = element->GetInput2(AdjointxEnum); _assert_(adjointx_input);
1515  Input2* adjointy_input = NULL;
1516  Input2* rheologyb_input = element->GetInput2(MaterialsRheologyBEnum); _assert_(rheologyb_input);
1517  if(domaintype!=Domain2DverticalEnum){
1518  vy_input = element->GetInput2(VyEnum); _assert_(vy_input);
1519  adjointy_input = element->GetInput2(AdjointyEnum); _assert_(adjointy_input);
1520  }
1521  /* Start looping on the number of gaussian points: */
1522  Gauss* gauss=element->NewGauss(4);
1523  for(int ig=gauss->begin();ig<gauss->end();ig++){
1524  gauss->GaussPoint(ig);
1525 
1526  thickness_input->GetInputValue(&thickness,gauss);
1527  vx_input->GetInputDerivativeValue(&dvx[0],xyz_list,gauss);
1528  adjointx_input->GetInputDerivativeValue(&dadjx[0],xyz_list,gauss);
1529  dim=2;
1530  if(domaintype!=Domain2DverticalEnum){
1531  adjointy_input->GetInputDerivativeValue(&dadjy[0],xyz_list, gauss);
1532  vy_input->GetInputDerivativeValue(&dvy[0],xyz_list,gauss);
1533  dim=3;
1534  }
1535 
1536  element->dViscositydBHO(&dmudB,dim,xyz_list,gauss,vx_input,vy_input);
1537 
1538  element->JacobianDeterminant(&Jdet,xyz_list,gauss);
1539  element->NodalFunctionsP1(basis,gauss);
1540 
1541  /*Build gradient vector (actually -dJ/dB): */
1542  for(int i=0;i<numvertices;i++){
1543  if(domaintype!=Domain2DverticalEnum){
1544  ge[i]+=-dmudB*thickness*(
1545  (2*dvx[0]+dvy[1])*2*dadjx[0]+(dvx[1]+dvy[0])*(dadjx[1]+dadjy[0])+(2*dvy[1]+dvx[0])*2*dadjy[1]
1546  )*Jdet*gauss->weight*basis[i];
1547  }
1548  else{
1549  ge[i]+=-dmudB*thickness*4*dvx[0]*dadjx[0]*Jdet*gauss->weight*basis[i];
1550  }
1551  _assert_(!xIsNan<IssmDouble>(ge[i]));
1552  }
1553  }
1554  gradient->SetValues(numvertices,vertexpidlist,ge,ADD_VAL);
1555 
1556  /*Clean up and return*/
1557  xDelete<IssmDouble>(xyz_list);
1558  xDelete<IssmDouble>(basis);
1559  xDelete<IssmDouble>(ge);
1560  xDelete<int>(vertexpidlist);
1561  delete gauss;
1562 }/*}}}*/

◆ GradientJBSSA()

void AdjointHorizAnalysis::GradientJBSSA ( Element *  element,
Vector< IssmDouble > *  gradient,
int  control_index 
)

Definition at line 1563 of file AdjointHorizAnalysis.cpp.

1563  {/*{{{*/
1564 
1565  /*Intermediaries*/
1566  int domaintype,dim;
1567  Element* basalelement;
1568 
1569  /*Get basal element*/
1570  element->FindParam(&domaintype,DomainTypeEnum);
1571  switch(domaintype){
1572  case Domain2DhorizontalEnum:
1573  basalelement = element;
1574  dim = 2;
1575  break;
1576  case Domain2DverticalEnum:
1577  if(!element->IsOnBase()) return;
1578  basalelement = element->SpawnBasalElement();
1579  dim = 1;
1580  break;
1581  case Domain3DEnum:
1582  if(!element->IsOnBase()) return;
1583  basalelement = element->SpawnBasalElement();
1584  dim = 2;
1585  break;
1586  default: _error_("mesh "<<EnumToStringx(domaintype)<<" not supported yet");
1587  }
1588 
1589  /*Intermediaries*/
1590  IssmDouble Jdet,weight;
1591  IssmDouble thickness,dmudB;
1592  IssmDouble dvx[3],dvy[3],dadjx[3],dadjy[3];
1593  IssmDouble *xyz_list= NULL;
1594 
1595  /*Fetch number of vertices for this finite element*/
1596  int numvertices = basalelement->GetNumberOfVertices();
1597 
1598  /*Initialize some vectors*/
1599  IssmDouble* basis = xNew<IssmDouble>(numvertices);
1600  IssmDouble* ge = xNewZeroInit<IssmDouble>(numvertices);
1601  int* vertexpidlist = xNew<int>(numvertices);
1602 
1603  /*Retrieve all inputs we will be needing: */
1604  basalelement->GetVerticesCoordinates(&xyz_list);
1605  basalelement->GradientIndexing(&vertexpidlist[0],control_index);
1606  Input2* thickness_input = basalelement->GetInput2(ThicknessEnum); _assert_(thickness_input);
1607  Input2* vx_input = basalelement->GetInput2(VxEnum); _assert_(vx_input);
1608  Input2* vy_input = basalelement->GetInput2(VyEnum); _assert_(vy_input);
1609  Input2* adjointx_input = basalelement->GetInput2(AdjointxEnum); _assert_(adjointx_input);
1610  Input2* adjointy_input = basalelement->GetInput2(AdjointyEnum); _assert_(adjointy_input);
1611  Input2* rheologyb_input = basalelement->GetInput2(MaterialsRheologyBEnum); _assert_(rheologyb_input);
1612 
1613  /* Start looping on the number of gaussian points: */
1614  Gauss* gauss=basalelement->NewGauss(4);
1615  for(int ig=gauss->begin();ig<gauss->end();ig++){
1616  gauss->GaussPoint(ig);
1617 
1618  thickness_input->GetInputValue(&thickness,gauss);
1619  vx_input->GetInputDerivativeValue(&dvx[0],xyz_list,gauss);
1620  vy_input->GetInputDerivativeValue(&dvy[0],xyz_list,gauss);
1621  adjointx_input->GetInputDerivativeValue(&dadjx[0],xyz_list,gauss);
1622  adjointy_input->GetInputDerivativeValue(&dadjy[0],xyz_list,gauss);
1623 
1624  basalelement->dViscositydBSSA(&dmudB,dim,xyz_list,gauss,vx_input,vy_input);
1625 
1626  basalelement->JacobianDeterminant(&Jdet,xyz_list,gauss);
1627  basalelement->NodalFunctionsP1(basis,gauss);
1628 
1629  /*Build gradient vector (actually -dJ/dB): */
1630  for(int i=0;i<numvertices;i++){
1631  ge[i]+=-dmudB*thickness*(
1632  (2*dvx[0]+dvy[1])*2*dadjx[0]+(dvx[1]+dvy[0])*(dadjx[1]+dadjy[0])+(2*dvy[1]+dvx[0])*2*dadjy[1]
1633  )*Jdet*gauss->weight*basis[i];
1634  _assert_(!xIsNan<IssmDouble>(ge[i]));
1635  }
1636  }
1637  gradient->SetValues(numvertices,vertexpidlist,ge,ADD_VAL);
1638 
1639  /*Clean up and return*/
1640  xDelete<IssmDouble>(xyz_list);
1641  xDelete<IssmDouble>(basis);
1642  xDelete<IssmDouble>(ge);
1643  xDelete<int>(vertexpidlist);
1644  delete gauss;
1645  if(domaintype!=Domain2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;};
1646 }/*}}}*/

◆ GradientJDragFS()

void AdjointHorizAnalysis::GradientJDragFS ( Element *  element,
Vector< IssmDouble > *  gradient,
int  control_index 
)

Definition at line 1729 of file AdjointHorizAnalysis.cpp.

1729  {/*{{{*/
1730 
1731  /*return if floating or not on bed (gradient is 0)*/
1732  if(element->IsFloating()) return;
1733  if(!element->IsOnBase()) return;
1734 
1735  /*Intermediaries*/
1736  int domaintype,dim;
1737  IssmDouble Jdet,weight;
1738  IssmDouble drag,dalpha2dk,normal[3];
1739  IssmDouble vx,vy,vz,lambda,mu,xi;
1740  IssmDouble *xyz_list_base= NULL;
1741 
1742  /*Fetch number of vertices for this finite element*/
1743  int numvertices = element->GetNumberOfVertices();
1744 
1745  /*Initialize some vectors*/
1746  IssmDouble* basis = xNew<IssmDouble>(numvertices);
1747  IssmDouble* ge = xNewZeroInit<IssmDouble>(numvertices);
1748  int* vertexpidlist = xNew<int>(numvertices);
1749 
1750  /* get domaintype */
1751  element->FindParam(&domaintype,DomainTypeEnum);
1752 
1753  /*Build friction element, needed later: */
1754  if(domaintype!=Domain2DverticalEnum) dim=3;
1755  else dim=2;
1756  Friction* friction=new Friction(element,dim);
1757 
1758  /*Retrieve all inputs we will be needing: */
1759  element->GetVerticesCoordinatesBase(&xyz_list_base);
1760  element->GradientIndexing(&vertexpidlist[0],control_index);
1761  Input2* vx_input = element->GetInput2(VxEnum); _assert_(vx_input);
1762  Input2* vy_input = element->GetInput2(VyEnum); _assert_(vy_input);
1763  Input2* adjointx_input = element->GetInput2(AdjointxEnum); _assert_(adjointx_input);
1764  Input2* adjointy_input = element->GetInput2(AdjointyEnum); _assert_(adjointy_input);
1765  Input2* vz_input = NULL;
1766  Input2* adjointz_input = NULL;
1767  if(domaintype!=Domain2DverticalEnum){
1768  vz_input = element->GetInput2(VzEnum); _assert_(vy_input);
1769  adjointz_input = element->GetInput2(AdjointzEnum); _assert_(adjointz_input);
1770  }
1771  Input2* dragcoeff_input = element->GetInput2(FrictionCoefficientEnum); _assert_(dragcoeff_input);
1772 
1773  /* Start looping on the number of gaussian points: */
1774  Gauss* gauss=element->NewGaussBase(4);
1775  for(int ig=gauss->begin();ig<gauss->end();ig++){
1776  gauss->GaussPoint(ig);
1777 
1778  adjointx_input->GetInputValue(&lambda, gauss);
1779  adjointy_input->GetInputValue(&mu, gauss);
1780  vx_input->GetInputValue(&vx,gauss);
1781  vy_input->GetInputValue(&vy,gauss);
1782  if(domaintype!=Domain2DverticalEnum){
1783  adjointz_input->GetInputValue(&xi ,gauss);
1784  vz_input->GetInputValue(&vz,gauss);
1785  }
1786  dragcoeff_input->GetInputValue(&drag, gauss);
1787 
1788  friction->GetAlphaComplement(&dalpha2dk,gauss);
1789  element->NormalBase(&normal[0],xyz_list_base);
1790 
1791  element->JacobianDeterminantBase(&Jdet,xyz_list_base,gauss);
1792  element->NodalFunctionsP1(basis,gauss);
1793 
1794  /*Build gradient vector (actually -dJ/dk): */
1795  if(domaintype!=Domain2DverticalEnum){
1796  for(int i=0;i<numvertices;i++){
1797  ge[i]+=(
1798  -lambda*(2*drag*dalpha2dk*(vx - vz*normal[0]*normal[2]))
1799  -mu *(2*drag*dalpha2dk*(vy - vz*normal[1]*normal[2]))
1800  -xi *(2*drag*dalpha2dk*(-vx*normal[0]*normal[2]-vy*normal[1]*normal[2]))
1801  )*Jdet*gauss->weight*basis[i];
1802  _assert_(!xIsNan<IssmDouble>(ge[i]));
1803  }
1804  }
1805  else{
1806  for(int i=0;i<numvertices;i++){
1807  ge[i]+=(
1808  -lambda*2*drag*dalpha2dk*vx
1809  -mu *2*drag*dalpha2dk*vy
1810  )*Jdet*gauss->weight*basis[i];
1811  _assert_(!xIsNan<IssmDouble>(ge[i]));
1812  }
1813  }
1814  }
1815  gradient->SetValues(numvertices,vertexpidlist,ge,ADD_VAL);
1816 
1817  /*Clean up and return*/
1818  xDelete<IssmDouble>(xyz_list_base);
1819  xDelete<IssmDouble>(basis);
1820  xDelete<IssmDouble>(ge);
1821  xDelete<int>(vertexpidlist);
1822  delete gauss;
1823  delete friction;
1824 }/*}}}*/

◆ GradientJDragGradient()

void AdjointHorizAnalysis::GradientJDragGradient ( Element *  element,
Vector< IssmDouble > *  gradient,
int  control_index 
)

Definition at line 1647 of file AdjointHorizAnalysis.cpp.

1647  {/*{{{*/
1648 
1649  /*return if floating (gradient is 0)*/
1650  if(element->IsFloating()) return;
1651 
1652  /*Intermediaries*/
1653  int domaintype,dim;
1654  Element* basalelement;
1655 
1656  /*Get basal element*/
1657  element->FindParam(&domaintype,DomainTypeEnum);
1658  switch(domaintype){
1659  case Domain2DhorizontalEnum:
1660  basalelement = element;
1661  dim = 2;
1662  break;
1663  case Domain2DverticalEnum:
1664  if(!element->IsOnBase()) return;
1665  basalelement = element->SpawnBasalElement();
1666  dim = 1;
1667  break;
1668  case Domain3DEnum:
1669  if(!element->IsOnBase()) return;
1670  basalelement = element->SpawnBasalElement();
1671  dim = 2;
1672  break;
1673  default: _error_("mesh "<<EnumToStringx(domaintype)<<" not supported yet");
1674  }
1675 
1676  /*Intermediaries*/
1677  IssmDouble Jdet,weight;
1678  IssmDouble dk[3];
1679  IssmDouble *xyz_list= NULL;
1680 
1681  /*Fetch number of vertices for this finite element*/
1682  int numvertices = basalelement->GetNumberOfVertices();
1683 
1684  /*Initialize some vectors*/
1685  IssmDouble* dbasis = xNew<IssmDouble>(2*numvertices);
1686  IssmDouble* ge = xNewZeroInit<IssmDouble>(numvertices);
1687  int* vertexpidlist = xNew<int>(numvertices);
1688 
1689  /*Retrieve all inputs we will be needing: */
1690  basalelement->GetVerticesCoordinates(&xyz_list);
1691  basalelement->GradientIndexing(&vertexpidlist[0],control_index);
1692  Input2* dragcoefficient_input = basalelement->GetInput2(FrictionCoefficientEnum); _assert_(dragcoefficient_input);
1693  DatasetInput2* weights_input = basalelement->GetDatasetInput2(InversionCostFunctionsCoefficientsEnum); _assert_(weights_input);
1694 
1695  /* Start looping on the number of gaussian points: */
1696  Gauss* gauss=basalelement->NewGauss(2);
1697  for(int ig=gauss->begin();ig<gauss->end();ig++){
1698  gauss->GaussPoint(ig);
1699 
1700  basalelement->JacobianDeterminant(&Jdet,xyz_list,gauss);
1701  basalelement->NodalFunctionsP1Derivatives(dbasis,xyz_list,gauss);
1702  weights_input->GetInputValue(&weight,gauss,DragCoefficientAbsGradientEnum);
1703 
1704  /*Build alpha_complement_list: */
1705  dragcoefficient_input->GetInputDerivativeValue(&dk[0],xyz_list,gauss);
1706 
1707  /*Build gradient vector (actually -dJ/ddrag): */
1708  for(int i=0;i<numvertices;i++){
1709  if(dim==2){
1710  ge[i]+=-weight*Jdet*gauss->weight*(dbasis[0*numvertices+i]*dk[0]+dbasis[1*numvertices+i]*dk[1]);
1711  }
1712  else{
1713  ge[i]+=-weight*Jdet*gauss->weight*dbasis[0*numvertices+i]*dk[0];
1714  }
1715  _assert_(!xIsNan<IssmDouble>(ge[i]));
1716  }
1717  }
1718  gradient->SetValues(numvertices,vertexpidlist,ge,ADD_VAL);
1719 
1720  /*Clean up and return*/
1721  xDelete<IssmDouble>(xyz_list);
1722  xDelete<IssmDouble>(dbasis);
1723  xDelete<IssmDouble>(ge);
1724  xDelete<int>(vertexpidlist);
1725  delete gauss;
1726  if(domaintype!=Domain2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;};
1727 
1728 }/*}}}*/

◆ GradientJDragL1L2()

void AdjointHorizAnalysis::GradientJDragL1L2 ( Element *  element,
Vector< IssmDouble > *  gradient,
int  control_index 
)

Definition at line 1825 of file AdjointHorizAnalysis.cpp.

1825  {/*{{{*/
1826 
1827  /*Same as SSA*/
1828  return this->GradientJDragSSA(element,gradient,control_index);
1829 }/*}}}*/

◆ GradientJDragHO()

void AdjointHorizAnalysis::GradientJDragHO ( Element *  element,
Vector< IssmDouble > *  gradient,
int  control_index 
)

Definition at line 1830 of file AdjointHorizAnalysis.cpp.

1830  {/*{{{*/
1831 
1832  /*return if floating or not on bed (gradient is 0)*/
1833  if(element->IsFloating()) return;
1834  if(!element->IsOnBase()) return;
1835 
1836  /*Intermediaries*/
1837  IssmDouble Jdet,weight;
1838  IssmDouble drag,dalpha2dk;
1839  IssmDouble vx,vy,lambda,mu;
1840  IssmDouble *xyz_list_base= NULL;
1841 
1842  int domaintype,dim;
1843  element->FindParam(&domaintype,DomainTypeEnum);
1844  /*Fetch number of vertices for this finite element*/
1845  int numvertices = element->GetNumberOfVertices();
1846 
1847  /*Initialize some vectors*/
1848  IssmDouble* basis = xNew<IssmDouble>(numvertices);
1849  IssmDouble* ge = xNewZeroInit<IssmDouble>(numvertices);
1850  int* vertexpidlist = xNew<int>(numvertices);
1851 
1852  /*Build friction element, needed later: */
1853  if(domaintype!=Domain2DverticalEnum) dim=3;
1854  else dim=2;
1855  Friction* friction=new Friction(element,dim);
1856 
1857  /*Retrieve all inputs we will be needing: */
1858  element->GetVerticesCoordinatesBase(&xyz_list_base);
1859  element->GradientIndexing(&vertexpidlist[0],control_index);
1860  Input2* vx_input = element->GetInput2(VxEnum); _assert_(vx_input);
1861  Input2* vy_input = NULL;
1862  Input2* adjointx_input = element->GetInput2(AdjointxEnum); _assert_(adjointx_input);
1863  Input2* adjointy_input = NULL;
1864  Input2* dragcoeff_input = element->GetInput2(FrictionCoefficientEnum); _assert_(dragcoeff_input);
1865  if(domaintype!=Domain2DverticalEnum){
1866  vy_input = element->GetInput2(VyEnum); _assert_(vy_input);
1867  adjointy_input = element->GetInput2(AdjointyEnum); _assert_(adjointy_input);
1868  }
1869  /* Start looping on the number of gaussian points: */
1870  Gauss* gauss=element->NewGaussBase(4);
1871  for(int ig=gauss->begin();ig<gauss->end();ig++){
1872  gauss->GaussPoint(ig);
1873 
1874  adjointx_input->GetInputValue(&lambda, gauss);
1875  vx_input->GetInputValue(&vx,gauss);
1876  if(domaintype!=Domain2DverticalEnum){
1877  adjointy_input->GetInputValue(&mu, gauss);
1878  vy_input->GetInputValue(&vy,gauss);
1879  }
1880  dragcoeff_input->GetInputValue(&drag, gauss);
1881 
1882  friction->GetAlphaComplement(&dalpha2dk,gauss);
1883 
1884  element->JacobianDeterminantBase(&Jdet,xyz_list_base,gauss);
1885  element->NodalFunctionsP1(basis,gauss);
1886 
1887  /*Build gradient vector (actually -dJ/dD): */
1888  for(int i=0;i<numvertices;i++){
1889  if(domaintype!=Domain2DverticalEnum) ge[i]+=-2.*drag*dalpha2dk*((lambda*vx+mu*vy))*Jdet*gauss->weight*basis[i];
1890  else ge[i]+=-2.*drag*dalpha2dk*(lambda*vx)*Jdet*gauss->weight*basis[i];
1891  _assert_(!xIsNan<IssmDouble>(ge[i]));
1892  }
1893  }
1894  gradient->SetValues(numvertices,vertexpidlist,ge,ADD_VAL);
1895 
1896  /*Clean up and return*/
1897  xDelete<IssmDouble>(xyz_list_base);
1898  xDelete<IssmDouble>(basis);
1899  xDelete<IssmDouble>(ge);
1900  xDelete<int>(vertexpidlist);
1901  delete gauss;
1902  delete friction;
1903 }/*}}}*/

◆ GradientJDragSSA()

void AdjointHorizAnalysis::GradientJDragSSA ( Element *  element,
Vector< IssmDouble > *  gradient,
int  control_index 
)

Definition at line 1904 of file AdjointHorizAnalysis.cpp.

1904  {/*{{{*/
1905 
1906  /*return if floating (gradient is 0)*/
1907  if(element->IsFloating()) return;
1908 
1909  /*Intermediaries*/
1910  int domaintype,dim;
1911  Element* basalelement;
1912 
1913  /*Get basal element*/
1914  element->FindParam(&domaintype,DomainTypeEnum);
1915  switch(domaintype){
1916  case Domain2DhorizontalEnum:
1917  basalelement = element;
1918  dim = 2;
1919  break;
1920  case Domain2DverticalEnum:
1921  if(!element->IsOnBase()) return;
1922  basalelement = element->SpawnBasalElement();
1923  dim = 1;
1924  break;
1925  case Domain3DEnum:
1926  if(!element->IsOnBase()) return;
1927  basalelement = element->SpawnBasalElement();
1928  dim = 2;
1929  break;
1930  default: _error_("mesh "<<EnumToStringx(domaintype)<<" not supported yet");
1931  }
1932 
1933  /*Intermediaries*/
1934  IssmDouble Jdet,weight;
1935  IssmDouble drag,dalpha2dk;
1936  IssmDouble vx,vy,lambda,mu;
1937  IssmDouble *xyz_list= NULL;
1938 
1939  /*Fetch number of vertices for this finite element*/
1940  int numvertices = basalelement->GetNumberOfVertices();
1941 
1942  /*Initialize some vectors*/
1943  IssmDouble* basis = xNew<IssmDouble>(numvertices);
1944  IssmDouble* ge = xNewZeroInit<IssmDouble>(numvertices);
1945  int* vertexpidlist = xNew<int>(numvertices);
1946 
1947  /*Build friction element, needed later: */
1948  Friction* friction=new Friction(basalelement,dim);
1949 
1950  /*Retrieve all inputs we will be needing: */
1951  basalelement->GetVerticesCoordinates(&xyz_list);
1952  basalelement->GradientIndexing(&vertexpidlist[0],control_index);
1953  Input2* vx_input = basalelement->GetInput2(VxEnum); _assert_(vx_input);
1954  Input2* vy_input = basalelement->GetInput2(VyEnum); _assert_(vy_input);
1955  Input2* adjointx_input = basalelement->GetInput2(AdjointxEnum); _assert_(adjointx_input);
1956  Input2* adjointy_input = basalelement->GetInput2(AdjointyEnum); _assert_(adjointy_input);
1957  Input2* dragcoeff_input = basalelement->GetInput2(FrictionCoefficientEnum); _assert_(dragcoeff_input);
1958 
1959  /* Start looping on the number of gaussian points: */
1960  Gauss* gauss=basalelement->NewGauss(4);
1961  for(int ig=gauss->begin();ig<gauss->end();ig++){
1962  gauss->GaussPoint(ig);
1963 
1964  adjointx_input->GetInputValue(&lambda, gauss);
1965  adjointy_input->GetInputValue(&mu, gauss);
1966  vx_input->GetInputValue(&vx,gauss);
1967  vy_input->GetInputValue(&vy,gauss);
1968  dragcoeff_input->GetInputValue(&drag, gauss);
1969 
1970  friction->GetAlphaComplement(&dalpha2dk,gauss);
1971 
1972  basalelement->JacobianDeterminant(&Jdet,xyz_list,gauss);
1973  basalelement->NodalFunctionsP1(basis,gauss);
1974 
1975  /*Build gradient vector (actually -dJ/dD): */
1976  for(int i=0;i<numvertices;i++){
1977  ge[i]+=-2.*drag*dalpha2dk*((lambda*vx+mu*vy))*Jdet*gauss->weight*basis[i];
1978  _assert_(!xIsNan<IssmDouble>(ge[i]));
1979  }
1980  }
1981  gradient->SetValues(numvertices,vertexpidlist,ge,ADD_VAL);
1982 
1983  /*Clean up and return*/
1984  xDelete<IssmDouble>(xyz_list);
1985  xDelete<IssmDouble>(basis);
1986  xDelete<IssmDouble>(ge);
1987  xDelete<int>(vertexpidlist);
1988  delete gauss;
1989  delete friction;
1990  if(domaintype!=Domain2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;};
1991 }/*}}}*/

◆ GradientJDragHydroFS()

void AdjointHorizAnalysis::GradientJDragHydroFS ( Element *  element,
Vector< IssmDouble > *  gradient,
int  control_index 
)

Definition at line 1993 of file AdjointHorizAnalysis.cpp.

1993  {/*{{{*/
1994 
1995  /*return if floating or not on bed (gradient is 0)*/
1996  if(element->IsFloating()) return;
1997  if(!element->IsOnBase()) return;
1998 
1999  /*Intermediaries*/
2000  int domaintype,dim;
2001  IssmDouble Jdet,weight;
2002  IssmDouble drag,dalpha2dk,normal[3];
2003  IssmDouble vx,vy,vz,lambda,mu,xi;
2004  IssmDouble *xyz_list_base= NULL;
2005 
2006  /*Fetch number of vertices for this finite element*/
2007  int numvertices = element->GetNumberOfVertices();
2008 
2009  /*Initialize some vectors*/
2010  IssmDouble* basis = xNew<IssmDouble>(numvertices);
2011  IssmDouble* ge = xNewZeroInit<IssmDouble>(numvertices);
2012  int* vertexpidlist = xNew<int>(numvertices);
2013 
2014  /* get domaintype */
2015  element->FindParam(&domaintype,DomainTypeEnum);
2016 
2017  /*Build friction element, needed later: */
2018  if(domaintype!=Domain2DverticalEnum) dim=3;
2019  else dim=2;
2020  Friction* friction=new Friction(element,dim);
2021 
2022  /*Retrieve all inputs we will be needing: */
2023  element->GetVerticesCoordinatesBase(&xyz_list_base);
2024  element->GradientIndexing(&vertexpidlist[0],control_index);
2025  Input2* vx_input = element->GetInput2(VxEnum); _assert_(vx_input);
2026  Input2* vy_input = element->GetInput2(VyEnum); _assert_(vy_input);
2027  Input2* adjointx_input = element->GetInput2(AdjointxEnum); _assert_(adjointx_input);
2028  Input2* adjointy_input = element->GetInput2(AdjointyEnum); _assert_(adjointy_input);
2029  Input2* vz_input = NULL;
2030  Input2* adjointz_input = NULL;
2031  if(domaintype!=Domain2DverticalEnum){
2032  vz_input = element->GetInput2(VzEnum); _assert_(vy_input);
2033  adjointz_input = element->GetInput2(AdjointzEnum); _assert_(adjointz_input);
2034  }
2035 
2036  /* Start looping on the number of gaussian points: */
2037  Gauss* gauss=element->NewGaussBase(4);
2038  for(int ig=gauss->begin();ig<gauss->end();ig++){
2039  gauss->GaussPoint(ig);
2040 
2041  adjointx_input->GetInputValue(&lambda, gauss);
2042  adjointy_input->GetInputValue(&mu, gauss);
2043  vx_input->GetInputValue(&vx,gauss);
2044  vy_input->GetInputValue(&vy,gauss);
2045  if(domaintype!=Domain2DverticalEnum){
2046  adjointz_input->GetInputValue(&xi ,gauss);
2047  vz_input->GetInputValue(&vz,gauss);
2048  }
2049 
2050  friction->GetAlphaComplement(&dalpha2dk,gauss);
2051  element->NormalBase(&normal[0],xyz_list_base);
2052 
2053  element->JacobianDeterminantBase(&Jdet,xyz_list_base,gauss);
2054  element->NodalFunctionsP1(basis,gauss);
2055 
2056  /*Build gradient vector (actually -dJ/dk): */
2057  if(domaintype!=Domain2DverticalEnum){
2058  for(int i=0;i<numvertices;i++){
2059  ge[i]+=(
2060  -lambda*(dalpha2dk*(vx - vz*normal[0]*normal[2]))
2061  -mu *(dalpha2dk*(vy - vz*normal[1]*normal[2]))
2062  -xi *(dalpha2dk*(-vx*normal[0]*normal[2]-vy*normal[1]*normal[2]))
2063  )*Jdet*gauss->weight*basis[i];
2064  _assert_(!xIsNan<IssmDouble>(ge[i]));
2065  }
2066  }
2067  else{
2068  for(int i=0;i<numvertices;i++){
2069  ge[i]+=(
2070  -lambda*dalpha2dk*vx
2071  -mu *dalpha2dk*vy
2072  )*Jdet*gauss->weight*basis[i];
2073  _assert_(!xIsNan<IssmDouble>(ge[i]));
2074  }
2075  }
2076  }
2077  gradient->SetValues(numvertices,vertexpidlist,ge,ADD_VAL);
2078 
2079  /*Clean up and return*/
2080  xDelete<IssmDouble>(xyz_list_base);
2081  xDelete<IssmDouble>(basis);
2082  xDelete<IssmDouble>(ge);
2083  xDelete<int>(vertexpidlist);
2084  delete gauss;
2085  delete friction;
2086 }/*}}}*/

◆ GradientJDragHydroL1L2()

void AdjointHorizAnalysis::GradientJDragHydroL1L2 ( Element *  element,
Vector< IssmDouble > *  gradient,
int  control_index 
)

Definition at line 2087 of file AdjointHorizAnalysis.cpp.

2087  {/*{{{*/
2088 
2089  /*Same as SSA*/
2090  return this->GradientJDragSSA(element,gradient,control_index);
2091 }/*}}}*/

◆ GradientJDragHydroHO()

void AdjointHorizAnalysis::GradientJDragHydroHO ( Element *  element,
Vector< IssmDouble > *  gradient,
int  control_index 
)

Definition at line 2092 of file AdjointHorizAnalysis.cpp.

2092  {/*{{{*/
2093 
2094  /*return if floating or not on bed (gradient is 0)*/
2095  if(element->IsFloating()) return;
2096  if(!element->IsOnBase()) return;
2097 
2098  /*Intermediaries*/
2099  IssmDouble Jdet,weight;
2100  IssmDouble drag,dalpha2dk;
2101  IssmDouble vx,vy,lambda,mu;
2102  IssmDouble *xyz_list_base= NULL;
2103 
2104  int domaintype,dim;
2105  element->FindParam(&domaintype,DomainTypeEnum);
2106  /*Fetch number of vertices for this finite element*/
2107  int numvertices = element->GetNumberOfVertices();
2108 
2109  /*Initialize some vectors*/
2110  IssmDouble* basis = xNew<IssmDouble>(numvertices);
2111  IssmDouble* ge = xNewZeroInit<IssmDouble>(numvertices);
2112  int* vertexpidlist = xNew<int>(numvertices);
2113 
2114  /*Build friction element, needed later: */
2115  if(domaintype!=Domain2DverticalEnum) dim=3;
2116  else dim=2;
2117  Friction* friction=new Friction(element,dim);
2118 
2119  /*Retrieve all inputs we will be needing: */
2120  element->GetVerticesCoordinatesBase(&xyz_list_base);
2121  element->GradientIndexing(&vertexpidlist[0],control_index);
2122  Input2* vx_input = element->GetInput2(VxEnum); _assert_(vx_input);
2123  Input2* vy_input = NULL;
2124  Input2* adjointx_input = element->GetInput2(AdjointxEnum); _assert_(adjointx_input);
2125  Input2* adjointy_input = NULL;
2126  if(domaintype!=Domain2DverticalEnum){
2127  vy_input = element->GetInput2(VyEnum); _assert_(vy_input);
2128  adjointy_input = element->GetInput2(AdjointyEnum); _assert_(adjointy_input);
2129  }
2130  /* Start looping on the number of gaussian points: */
2131  Gauss* gauss=element->NewGaussBase(4);
2132  for(int ig=gauss->begin();ig<gauss->end();ig++){
2133  gauss->GaussPoint(ig);
2134 
2135  adjointx_input->GetInputValue(&lambda, gauss);
2136  vx_input->GetInputValue(&vx,gauss);
2137  if(domaintype!=Domain2DverticalEnum){
2138  adjointy_input->GetInputValue(&mu, gauss);
2139  vy_input->GetInputValue(&vy,gauss);
2140  }
2141 
2142  friction->GetAlphaComplement(&dalpha2dk,gauss);
2143 
2144  element->JacobianDeterminantBase(&Jdet,xyz_list_base,gauss);
2145  element->NodalFunctionsP1(basis,gauss);
2146 
2147  /*Build gradient vector (actually -dJ/dD): */
2148  for(int i=0;i<numvertices;i++){
2149  if(domaintype!=Domain2DverticalEnum) ge[i]+=-dalpha2dk*((lambda*vx+mu*vy))*Jdet*gauss->weight*basis[i];
2150  else ge[i]+=-dalpha2dk*(lambda*vx)*Jdet*gauss->weight*basis[i];
2151  _assert_(!xIsNan<IssmDouble>(ge[i]));
2152  }
2153  }
2154  gradient->SetValues(numvertices,vertexpidlist,ge,ADD_VAL);
2155 
2156  /*Clean up and return*/
2157  xDelete<IssmDouble>(xyz_list_base);
2158  xDelete<IssmDouble>(basis);
2159  xDelete<IssmDouble>(ge);
2160  xDelete<int>(vertexpidlist);
2161  delete gauss;
2162  delete friction;
2163 }/*}}}*/

◆ GradientJDragHydroSSA()

void AdjointHorizAnalysis::GradientJDragHydroSSA ( Element *  element,
Vector< IssmDouble > *  gradient,
int  control_index 
)

Definition at line 2165 of file AdjointHorizAnalysis.cpp.

2165  {/*{{{*/
2166 
2167  /*return if floating (gradient is 0)*/
2168  if(element->IsFloating()) return;
2169 
2170  /*Intermediaries*/
2171  int domaintype,dim;
2172 
2173  Element* basalelement;
2174 
2175  /*Get basal element*/
2176  element->FindParam(&domaintype,DomainTypeEnum);
2177  switch(domaintype){
2178  case Domain2DhorizontalEnum:
2179  basalelement = element;
2180  dim = 2;
2181  break;
2182  case Domain2DverticalEnum:
2183  if(!element->IsOnBase()) return;
2184  basalelement = element->SpawnBasalElement();
2185  dim = 1;
2186  break;
2187  case Domain3DEnum:
2188  if(!element->IsOnBase()) return;
2189  basalelement = element->SpawnBasalElement();
2190  dim = 2;
2191  break;
2192  default: _error_("mesh "<<EnumToStringx(domaintype)<<" not supported yet");
2193  }
2194 
2195  /*Intermediaries*/
2196  IssmDouble Jdet,weight;
2197  IssmDouble dalpha2dk;
2198  IssmDouble vx,vy,lambda,mu;
2199  IssmDouble *xyz_list= NULL;
2200 
2201  /*Fetch number of vertices for this finite element*/
2202  int numvertices = basalelement->GetNumberOfVertices();
2203 
2204  /*Initialize some vectors*/
2205  IssmDouble* basis = xNew<IssmDouble>(numvertices);
2206  IssmDouble* ge = xNewZeroInit<IssmDouble>(numvertices);
2207  int* vertexpidlist = xNew<int>(numvertices);
2208 
2209  /*Build friction element, needed later: */
2210  Friction* friction=new Friction(basalelement,dim);
2211 
2212  /*Retrieve all inputs we will be needing: */
2213  basalelement->GetVerticesCoordinates(&xyz_list);
2214  basalelement->GradientIndexing(&vertexpidlist[0],control_index);
2215  Input2* vx_input = basalelement->GetInput2(VxEnum); _assert_(vx_input);
2216  Input2* vy_input = basalelement->GetInput2(VyEnum); _assert_(vy_input);
2217  Input2* adjointx_input = basalelement->GetInput2(AdjointxEnum); _assert_(adjointx_input);
2218  Input2* adjointy_input = basalelement->GetInput2(AdjointyEnum); _assert_(adjointy_input);
2219 
2220  IssmDouble q_exp;
2221  IssmDouble C_param;
2222  IssmDouble As;
2223  IssmDouble Neff;
2224  IssmDouble n;
2225  IssmDouble alpha;
2226  IssmDouble Chi,Gamma;
2227  IssmDouble vz,vmag;
2228  IssmDouble Uder;
2229 
2230  /*Recover parameters: */
2231  Input2* qinput = basalelement->GetInput2(FrictionQEnum);
2232  Input2* cinput = basalelement->GetInput2(FrictionCEnum);
2233  Input2* Asinput = basalelement->GetInput2(FrictionAsEnum);
2234  Input2* nInput =basalelement->GetInput2(MaterialsRheologyNEnum);
2235  Input2* Ninput = basalelement->GetInput2(FrictionEffectivePressureEnum);
2236  /* Start looping on the number of gaussian points: */
2237  Gauss* gauss=basalelement->NewGauss(4);
2238  for(int ig=gauss->begin();ig<gauss->end();ig++){
2239  gauss->GaussPoint(ig);
2240 
2241  adjointx_input->GetInputValue(&lambda, gauss);
2242  adjointy_input->GetInputValue(&mu, gauss);
2243  vx_input->GetInputValue(&vx,gauss);
2244  vy_input->GetInputValue(&vy,gauss);
2245 
2246  friction->GetAlphaComplement(&dalpha2dk,gauss);
2247 
2248  basalelement->JacobianDeterminant(&Jdet,xyz_list,gauss);
2249  basalelement->NodalFunctionsP1(basis,gauss);
2250 
2251  /*Dealing with dalpha/du*/
2252  qinput->GetInputValue(&q_exp,gauss);
2253  cinput->GetInputValue(&C_param,gauss);
2254  Asinput->GetInputValue(&As,gauss);
2255  Ninput->GetInputValue(&Neff,gauss);
2256  nInput->GetInputValue(&n,gauss);
2257 
2258  if (q_exp==1){
2259  alpha=1;
2260  }
2261  else{
2262  alpha=(pow(q_exp-1,q_exp-1))/pow(q_exp,q_exp);
2263  }
2264 
2265  vmag = sqrt(vx*vx + vy*vy);
2266  Chi = vmag/(pow(C_param,n)*pow(Neff,n)*As);
2267  Gamma = (Chi/(1.+alpha*pow(Chi,q_exp)));
2268 
2269  Uder =Neff*C_param/(vmag*vmag*n) *
2270  (Gamma-alpha*q_exp*pow(Chi,q_exp-1.)*Gamma*Gamma* pow(Gamma,(1.-n)/n) -
2271  n* pow(Gamma,1./n));
2272 
2273  /*Build gradient vector (actually -dJ/dD): */
2274  for(int i=0;i<numvertices;i++){
2275  ge[i]+=-dalpha2dk*((lambda*vx+mu*vy))*Jdet*gauss->weight*basis[i];
2276  _assert_(!xIsNan<IssmDouble>(ge[i]));
2277  }
2278  }
2279  gradient->SetValues(numvertices,vertexpidlist,ge,ADD_VAL);
2280 
2281  /*Clean up and return*/
2282  xDelete<IssmDouble>(xyz_list);
2283  xDelete<IssmDouble>(basis);
2284  xDelete<IssmDouble>(ge);
2285  xDelete<int>(vertexpidlist);
2286  delete gauss;
2287  delete friction;
2288  if(domaintype!=Domain2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;};
2289 }/*}}}*/

◆ GradientJDSSA()

void AdjointHorizAnalysis::GradientJDSSA ( Element *  element,
Vector< IssmDouble > *  gradient,
int  control_index 
)

Definition at line 2291 of file AdjointHorizAnalysis.cpp.

2291  {/*{{{*/
2292 
2293  /*Intermediaries*/
2294  int domaintype,dim;
2295  Element* basalelement;
2296 
2297  /*Get basal element*/
2298  element->FindParam(&domaintype,DomainTypeEnum);
2299  switch(domaintype){
2300  case Domain2DhorizontalEnum:
2301  basalelement = element;
2302  dim = 2;
2303  break;
2304  case Domain2DverticalEnum:
2305  if(!element->IsOnBase()) return;
2306  basalelement = element->SpawnBasalElement();
2307  dim = 1;
2308  break;
2309  case Domain3DEnum:
2310  if(!element->IsOnBase()) return;
2311  basalelement = element->SpawnBasalElement();
2312  dim = 2;
2313  break;
2314  default: _error_("mesh "<<EnumToStringx(domaintype)<<" not supported yet");
2315  }
2316 
2317  /*Intermediaries*/
2318  IssmDouble Jdet,weight;
2319  IssmDouble thickness,dmudD;
2320  IssmDouble dvx[3],dvy[3],dadjx[3],dadjy[3];
2321  IssmDouble *xyz_list= NULL;
2322 
2323  /*Fetch number of vertices for this finite element*/
2324  int numvertices = basalelement->GetNumberOfVertices();
2325 
2326  /*Initialize some vectors*/
2327  IssmDouble* basis = xNew<IssmDouble>(numvertices);
2328  IssmDouble* ge = xNewZeroInit<IssmDouble>(numvertices);
2329  int* vertexpidlist = xNew<int>(numvertices);
2330 
2331  /*Retrieve all inputs we will be needing: */
2332  basalelement->GetVerticesCoordinates(&xyz_list);
2333  basalelement->GradientIndexing(&vertexpidlist[0],control_index);
2334  Input2* thickness_input = basalelement->GetInput2(ThicknessEnum); _assert_(thickness_input);
2335  Input2* vx_input = basalelement->GetInput2(VxEnum); _assert_(vx_input);
2336  Input2* vy_input = basalelement->GetInput2(VyEnum); _assert_(vy_input);
2337  Input2* adjointx_input = basalelement->GetInput2(AdjointxEnum); _assert_(adjointx_input);
2338  Input2* adjointy_input = basalelement->GetInput2(AdjointyEnum); _assert_(adjointy_input);
2339  Input2* rheologyb_input = basalelement->GetInput2(MaterialsRheologyBbarEnum); _assert_(rheologyb_input);
2340 
2341  /* Start looping on the number of gaussian points: */
2342  Gauss* gauss=basalelement->NewGauss(4);
2343  for(int ig=gauss->begin();ig<gauss->end();ig++){
2344  gauss->GaussPoint(ig);
2345 
2346  thickness_input->GetInputValue(&thickness,gauss);
2347  vx_input->GetInputDerivativeValue(&dvx[0],xyz_list,gauss);
2348  vy_input->GetInputDerivativeValue(&dvy[0],xyz_list,gauss);
2349  adjointx_input->GetInputDerivativeValue(&dadjx[0],xyz_list,gauss);
2350  adjointy_input->GetInputDerivativeValue(&dadjy[0],xyz_list,gauss);
2351 
2352  basalelement->dViscositydDSSA(&dmudD,dim,xyz_list,gauss,vx_input,vy_input);
2353 
2354  basalelement->JacobianDeterminant(&Jdet,xyz_list,gauss);
2355  basalelement->NodalFunctionsP1(basis,gauss);
2356 
2357  for(int i=0;i<numvertices;i++){
2358  ge[i]+=-dmudD*thickness*(
2359  (2*dvx[0]+dvy[1])*2*dadjx[0]+(dvx[1]+dvy[0])*(dadjx[1]+dadjy[0])+(2*dvy[1]+dvx[0])*2*dadjy[1]
2360  )*Jdet*gauss->weight*basis[i];
2361  _assert_(!xIsNan<IssmDouble>(ge[i]));
2362  }
2363  }
2364  gradient->SetValues(numvertices,vertexpidlist,ge,ADD_VAL);
2365 
2366  /*Clean up and return*/
2367  xDelete<IssmDouble>(xyz_list);
2368  xDelete<IssmDouble>(basis);
2369  xDelete<IssmDouble>(ge);
2370  xDelete<int>(vertexpidlist);
2371  delete gauss;
2372  if(domaintype!=Domain2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;};
2373 }/*}}}*/

◆ InputUpdateFromSolution()

void AdjointHorizAnalysis::InputUpdateFromSolution ( IssmDouble *  solution,
Element *  element 
)
virtual

Implements Analysis.

Definition at line 2374 of file AdjointHorizAnalysis.cpp.

2374  {/*{{{*/
2375  int approximation;
2376  element->GetInputValue(&approximation,ApproximationEnum);
2377  if(approximation==FSApproximationEnum || approximation==NoneApproximationEnum){
2378  InputUpdateFromSolutionFS(solution,element);
2379  }
2380  else{
2381  InputUpdateFromSolutionHoriz(solution,element);
2382  }
2383 }/*}}}*/

◆ InputUpdateFromSolutionFS()

void AdjointHorizAnalysis::InputUpdateFromSolutionFS ( IssmDouble *  solution,
Element *  element 
)

Definition at line 2384 of file AdjointHorizAnalysis.cpp.

2384  {/*{{{*/
2385  int i,fe_FS;
2386  int* vdoflist=NULL;
2387  int* pdoflist=NULL;
2388  IssmDouble FSreconditioning;
2389 
2390  int domaintype,dim;
2391  element->FindParam(&domaintype,DomainTypeEnum);
2392  switch(domaintype){
2393  case Domain2DhorizontalEnum:
2394  dim = 3;
2395  break;
2396  case Domain2DverticalEnum:
2397  dim = 2;
2398  break;
2399  case Domain3DEnum:
2400  dim = 3;
2401  break;
2402  default: _error_("mesh "<<EnumToStringx(domaintype)<<" not supported yet");
2403  }
2404 
2405  /*Fetch number of nodes and dof for this finite element*/
2406  int vnumnodes = element->NumberofNodesVelocity();
2407  int pnumnodes = element->NumberofNodesPressure();
2408  int vnumdof = vnumnodes*dim;
2409  int pnumdof = pnumnodes*1;
2410 
2411  /*Initialize values*/
2412  IssmDouble* values = xNew<IssmDouble>(vnumdof+pnumdof);
2413  IssmDouble* lambdax = xNew<IssmDouble>(vnumnodes);
2414  IssmDouble* lambday = xNew<IssmDouble>(vnumnodes);
2415  IssmDouble* lambdaz = xNew<IssmDouble>(vnumnodes);
2416  IssmDouble* lambdap = xNew<IssmDouble>(pnumnodes);
2417 
2418  int* cs_list = xNew<int>(vnumnodes+pnumnodes);
2419  if(dim==2) for(i=0;i<vnumnodes;i++) cs_list[i] = XYEnum;
2420  else for(i=0;i<vnumnodes;i++) cs_list[i] = XYZEnum;
2421  for(i=0;i<pnumnodes;i++) cs_list[vnumnodes+i] = PressureEnum;
2422 
2423  /*Get dof list: */
2424  element->GetDofListLocalVelocity(&vdoflist,GsetEnum);
2425  element->GetDofListLocalPressure(&pdoflist,GsetEnum);
2426 
2427  /*Use the dof list to index into the solution vector: */
2428  for(i=0;i<vnumdof;i++) values[i] =solution[vdoflist[i]];
2429  for(i=0;i<pnumdof;i++) values[vnumdof+i]=solution[pdoflist[i]];
2430 
2431  /*Transform solution in Cartesian Space*/
2432  element->TransformSolutionCoord(values,cs_list);
2433 
2434  /*fill in all arrays: */
2435  for(i=0;i<vnumnodes;i++){
2436  lambdax[i] = values[i*dim+0];
2437  if(xIsNan<IssmDouble>(lambdax[i])) _error_("NaN found in solution vector");
2438  if(xIsInf<IssmDouble>(lambdax[i])) _error_("Inf found in solution vector");
2439  lambday[i] = values[i*dim+1];
2440  if(xIsNan<IssmDouble>(lambday[i])) _error_("NaN found in solution vector");
2441  if(xIsInf<IssmDouble>(lambday[i])) _error_("Inf found in solution vector");
2442  if(dim==3){
2443  lambdaz[i] = values[i*dim+2];
2444  if(xIsNan<IssmDouble>(lambdaz[i])) _error_("NaN found in solution vector");
2445  if(xIsInf<IssmDouble>(lambdaz[i])) _error_("Inf found in solution vector");
2446  }
2447  }
2448  for(i=0;i<pnumnodes;i++){
2449  lambdap[i] = values[vnumdof+i];
2450  if(xIsNan<IssmDouble>(lambdap[i])) _error_("NaN found in solution vector");
2451  if(xIsInf<IssmDouble>(lambdap[i])) _error_("Inf found in solution vector");
2452  }
2453 
2454  /*Recondition pressure and compute vel: */
2455  element->FindParam(&FSreconditioning,StressbalanceFSreconditioningEnum);
2456  for(i=0;i<pnumnodes;i++) lambdap[i]=lambdap[i]*FSreconditioning;
2457 
2458  /*Add vx and vy as inputs to the tria element: */
2459  element->AddInput2(AdjointxEnum,lambdax,element->VelocityInterpolation());
2460  element->AddInput2(AdjointyEnum,lambday,element->VelocityInterpolation());
2461  if(domaintype!=Domain2DverticalEnum) element->AddInput2(AdjointzEnum,lambdaz,element->VelocityInterpolation());
2462 
2463  element->FindParam(&fe_FS,FlowequationFeFSEnum);
2464  if(fe_FS!=LATaylorHoodEnum && fe_FS!=LACrouzeixRaviartEnum)
2465  element->AddInput2(AdjointpEnum,lambdap,element->PressureInterpolation());
2466 
2467  /*Free ressources:*/
2468  xDelete<int>(vdoflist);
2469  xDelete<int>(pdoflist);
2470  xDelete<int>(cs_list);
2471  xDelete<IssmDouble>(lambdap);
2472  xDelete<IssmDouble>(lambdaz);
2473  xDelete<IssmDouble>(lambday);
2474  xDelete<IssmDouble>(lambdax);
2475  xDelete<IssmDouble>(values);
2476 }/*}}}*/

◆ InputUpdateFromSolutionHoriz()

void AdjointHorizAnalysis::InputUpdateFromSolutionHoriz ( IssmDouble *  solution,
Element *  element 
)

Definition at line 2477 of file AdjointHorizAnalysis.cpp.

2477  {/*{{{*/
2478  int i;
2479  int* doflist=NULL;
2480 
2481  int domaintype;
2482  element->FindParam(&domaintype,DomainTypeEnum);
2483 
2484  /*Fetch number of nodes and dof for this finite element*/
2485  int numnodes = element->GetNumberOfNodes();
2486  int numdof;
2487  if(domaintype!=Domain2DverticalEnum) numdof = numnodes*2;
2488  else numdof = numnodes*1;
2489  /*Fetch dof list and allocate solution vectors*/
2490  element->GetDofListLocal(&doflist,NoneApproximationEnum,GsetEnum);
2491  IssmDouble* values = xNew<IssmDouble>(numdof);
2492  IssmDouble* lambdax = xNew<IssmDouble>(numnodes);
2493  IssmDouble* lambday = xNew<IssmDouble>(numnodes);
2494 
2495  /*Use the dof list to index into the solution vector: */
2496  for(i=0;i<numdof;i++) values[i]=solution[doflist[i]];
2497 
2498  /*Transform solution in Cartesian Space*/
2499  if(domaintype!=Domain2DverticalEnum) element->TransformSolutionCoord(&values[0],XYEnum);
2500 
2501  /*Ok, we have vx and vy in values, fill in vx and vy arrays: */
2502  for(i=0;i<numnodes;i++){
2503  if(domaintype!=Domain2DverticalEnum){
2504  lambdax[i]=values[i*2+0];
2505  lambday[i]=values[i*2+1];
2506  }
2507  else {lambdax[i]=values[i];lambday[i]=0;}
2508  /*Check solution*/
2509  if(xIsNan<IssmDouble>(lambdax[i])) _error_("NaN found in solution vector");
2510  if(xIsInf<IssmDouble>(lambdax[i])) _error_("Inf found in solution vector");
2511  if(domaintype!=Domain2DverticalEnum && xIsNan<IssmDouble>(lambday[i])) _error_("NaN found in solution vector");
2512  if(domaintype!=Domain2DverticalEnum && xIsInf<IssmDouble>(lambday[i])) _error_("Inf found in solution vector");
2513  }
2514 
2515  /*Add vx and vy as inputs to the tria element: */
2516  element->AddInput2(AdjointxEnum,lambdax,element->GetElementType());
2517  element->AddInput2(AdjointyEnum,lambday,element->GetElementType());
2518 
2519  /*Free ressources:*/
2520  xDelete<IssmDouble>(values);
2521  xDelete<IssmDouble>(lambdax);
2522  xDelete<IssmDouble>(lambday);
2523  xDelete<int>(doflist);
2524 }/*}}}*/

◆ UpdateConstraints()

void AdjointHorizAnalysis::UpdateConstraints ( FemModel *  femmodel)
virtual

Implements Analysis.

Definition at line 2525 of file AdjointHorizAnalysis.cpp.

2525  {/*{{{*/
2526  /*Default, do nothing*/
2527  return;
2528 }/*}}}*/
The documentation for this class was generated from the following files:
FrictionCoefficientEnum
@ FrictionCoefficientEnum
Definition: EnumDefinitions.h:574
Element::TransformStiffnessMatrixCoord
void TransformStiffnessMatrixCoord(ElementMatrix *Ke, int cs_enum)
Definition: Element.cpp:4473
StressbalanceAnalysis
Definition: StressbalanceAnalysis.h:11
Element::GetElementType
virtual int GetElementType(void)=0
AdjointHorizAnalysis::GradientJBbarFS
void GradientJBbarFS(Element *element, Vector< IssmDouble > *gradient, int control_index)
Definition: AdjointHorizAnalysis.cpp:1180
StressbalanceAnalysis::CreateKMatrix
ElementMatrix * CreateKMatrix(Element *element)
Definition: StressbalanceAnalysis.cpp:1115
_assert_
#define _assert_(ignore)
Definition: exceptions.h:37
AdjointHorizAnalysis::GradientJDragGradient
void GradientJDragGradient(Element *element, Vector< IssmDouble > *gradient, int control_index)
Definition: AdjointHorizAnalysis.cpp:1647
IssmDouble
double IssmDouble
Definition: types.h:37
Element::TransformSolutionCoord
void TransformSolutionCoord(IssmDouble *solution, int cs_enum)
Definition: Element.cpp:4397
HOApproximationEnum
@ HOApproximationEnum
Definition: EnumDefinitions.h:1095
Element::IsOnBase
bool IsOnBase()
Definition: Element.cpp:1984
DatasetInput2
Definition: DatasetInput2.h:14
Element::GetNumberOfNodes
virtual int GetNumberOfNodes(void)=0
AdjointHorizAnalysis::CreateKMatrixHO
ElementMatrix * CreateKMatrixHO(Element *element)
Definition: AdjointHorizAnalysis.cpp:149
FrictionEffectivePressureEnum
@ FrictionEffectivePressureEnum
Definition: EnumDefinitions.h:576
Element::FindParam
void FindParam(bool *pvalue, int paramenum)
Definition: Element.cpp:933
AdjointHorizAnalysis::InputUpdateFromSolutionHoriz
void InputUpdateFromSolutionHoriz(IssmDouble *solution, Element *element)
Definition: AdjointHorizAnalysis.cpp:2477
InversionVyObsEnum
@ InversionVyObsEnum
Definition: EnumDefinitions.h:634
Element::dViscositydDSSA
void dViscositydDSSA(IssmDouble *pdmudB, int dim, IssmDouble *xyz_list, Gauss *gauss, Input2 *vx_input, Input2 *vy_input)
Definition: Element.cpp:874
Element::NewGaussBase
virtual Gauss * NewGaussBase(int order)=0
AdjointHorizAnalysis::GradientJBHO
void GradientJBHO(Element *element, Vector< IssmDouble > *gradient, int control_index)
Definition: AdjointHorizAnalysis.cpp:1487
InversionCostFunctionsCoefficientsEnum
@ InversionCostFunctionsCoefficientsEnum
Definition: EnumDefinitions.h:629
Friction::GetAlphaComplement
void GetAlphaComplement(IssmDouble *alpha_complement, Gauss *gauss)
Definition: Friction.cpp:42
StressbalanceFSreconditioningEnum
@ StressbalanceFSreconditioningEnum
Definition: EnumDefinitions.h:405
ADD_VAL
@ ADD_VAL
Definition: toolkitsenums.h:14
AdjointHorizAnalysis::InputUpdateFromSolutionFS
void InputUpdateFromSolutionFS(IssmDouble *solution, Element *element)
Definition: AdjointHorizAnalysis.cpp:2384
Element::PressureInterpolation
virtual int PressureInterpolation()=0
AdjointHorizAnalysis::CreateKMatrixL1L2
ElementMatrix * CreateKMatrixL1L2(Element *element)
Definition: AdjointHorizAnalysis.cpp:217
Element::SpawnBasalElement
virtual Element * SpawnBasalElement(void)=0
AdjointHorizAnalysis::GradientJDragHO
void GradientJDragHO(Element *element, Vector< IssmDouble > *gradient, int control_index)
Definition: AdjointHorizAnalysis.cpp:1830
Element::dViscositydBHO
void dViscositydBHO(IssmDouble *pdmudB, int dim, IssmDouble *xyz_list, Gauss *gauss, Input2 *vx_input, Input2 *vy_input)
Definition: Element.cpp:800
AdjointHorizAnalysis::GradientJDragL1L2
void GradientJDragL1L2(Element *element, Vector< IssmDouble > *gradient, int control_index)
Definition: AdjointHorizAnalysis.cpp:1825
Material::ViscosityHODerivativeEpsSquare
virtual void ViscosityHODerivativeEpsSquare(IssmDouble *pmu_prime, IssmDouble *epsilon, Gauss *gauss)=0
Element::TransformLoadVectorCoord
void TransformLoadVectorCoord(ElementVector *pe, int cs_enum)
Definition: Element.cpp:4346
PressureEnum
@ PressureEnum
Definition: EnumDefinitions.h:664
XYZEnum
@ XYZEnum
Definition: EnumDefinitions.h:1331
AdjointHorizAnalysis::GradientJDragHydroFS
void GradientJDragHydroFS(Element *element, Vector< IssmDouble > *gradient, int control_index)
Definition: AdjointHorizAnalysis.cpp:1993
ElementVector::values
IssmDouble * values
Definition: ElementVector.h:24
SurfaceRelVelMisfitEnum
@ SurfaceRelVelMisfitEnum
Definition: EnumDefinitions.h:828
AdjointHorizAnalysis::GradientJDragFS
void GradientJDragFS(Element *element, Vector< IssmDouble > *gradient, int control_index)
Definition: AdjointHorizAnalysis.cpp:1729
FrictionQEnum
@ FrictionQEnum
Definition: EnumDefinitions.h:580
Element::JacobianDeterminantBase
virtual void JacobianDeterminantBase(IssmDouble *Jdet, IssmDouble *xyz_list_base, Gauss *gauss)=0
MaterialsRheologyNEnum
@ MaterialsRheologyNEnum
Definition: EnumDefinitions.h:651
Element::JacobianDeterminantTop
virtual void JacobianDeterminantTop(IssmDouble *Jdet, IssmDouble *xyz_list_base, Gauss *gauss)=0
Element::IsFloating
bool IsFloating()
Definition: Element.cpp:1987
Element::NormalBase
virtual void NormalBase(IssmDouble *normal, IssmDouble *xyz_list)=0
Material::ViscosityFSDerivativeEpsSquare
virtual void ViscosityFSDerivativeEpsSquare(IssmDouble *pmu_prime, IssmDouble *epsilon, Gauss *gauss)=0
DatasetInput2::GetInputValue
void GetInputValue(IssmDouble *pvalue, Gauss *gauss, int index)
Definition: DatasetInput2.cpp:199
Element::NodalFunctionsP1Derivatives
virtual void NodalFunctionsP1Derivatives(IssmDouble *dbasis, IssmDouble *xyz_list, Gauss *gauss)=0
DragCoefficientAbsGradientEnum
@ DragCoefficientAbsGradientEnum
Definition: EnumDefinitions.h:537
FSApproximationEnum
@ FSApproximationEnum
Definition: EnumDefinitions.h:1060
VyEnum
@ VyEnum
Definition: EnumDefinitions.h:850
FlowequationFeFSEnum
@ FlowequationFeFSEnum
Definition: EnumDefinitions.h:137
FrictionAsEnum
@ FrictionAsEnum
Definition: EnumDefinitions.h:571
L1L2ApproximationEnum
@ L1L2ApproximationEnum
Definition: EnumDefinitions.h:1135
SurfaceLogVxVyMisfitEnum
@ SurfaceLogVxVyMisfitEnum
Definition: EnumDefinitions.h:826
Input2::GetInputDerivativeValue
virtual void GetInputDerivativeValue(IssmDouble *derivativevalues, IssmDouble *xyz_list, Gauss *gauss)
Definition: Input2.h:37
AdjointpEnum
@ AdjointpEnum
Definition: EnumDefinitions.h:465
MaterialsRheologyBbarEnum
@ MaterialsRheologyBbarEnum
Definition: EnumDefinitions.h:644
Element::GetDofListLocalPressure
void GetDofListLocalPressure(int **pdoflist, int setenum)
Definition: Element.cpp:1054
SurfaceAbsVelMisfitEnum
@ SurfaceAbsVelMisfitEnum
Definition: EnumDefinitions.h:818
Element::GetInput2
virtual Input2 * GetInput2(int inputenum)=0
Element::DeleteMaterials
void DeleteMaterials(void)
Definition: Element.cpp:429
Element
Definition: Element.h:41
Element::NodalFunctions
virtual void NodalFunctions(IssmDouble *basis, Gauss *gauss)=0
AdjointHorizAnalysis::GradientJDragHydroL1L2
void GradientJDragHydroL1L2(Element *element, Vector< IssmDouble > *gradient, int control_index)
Definition: AdjointHorizAnalysis.cpp:2087
Element::GetDatasetInput2
virtual DatasetInput2 * GetDatasetInput2(int inputenum)
Definition: Element.h:250
Domain2DhorizontalEnum
@ Domain2DhorizontalEnum
Definition: EnumDefinitions.h:534
RheologyBInitialguessEnum
@ RheologyBInitialguessEnum
Definition: EnumDefinitions.h:672
InversionIncompleteAdjointEnum
@ InversionIncompleteAdjointEnum
Definition: EnumDefinitions.h:217
Element::NodalFunctionsP1
virtual void NodalFunctionsP1(IssmDouble *basis, Gauss *gauss)=0
Element::NewElementVector
ElementVector * NewElementVector(int approximation_enum=NoneApproximationEnum)
Definition: Element.cpp:2505
LATaylorHoodEnum
@ LATaylorHoodEnum
Definition: EnumDefinitions.h:1139
XYEnum
@ XYEnum
Definition: EnumDefinitions.h:1330
NoneApproximationEnum
@ NoneApproximationEnum
Definition: EnumDefinitions.h:1201
Vector::SetValues
void SetValues(int ssize, int *list, doubletype *values, InsMode mode)
Definition: Vector.h:153
DomainTypeEnum
@ DomainTypeEnum
Definition: EnumDefinitions.h:124
Element::AddInput2
virtual void AddInput2(int input_enum, IssmDouble *values, int interpolation_enum)
Definition: Element.h:216
AdjointyEnum
@ AdjointyEnum
Definition: EnumDefinitions.h:467
Element::NewGauss
virtual Gauss * NewGauss(void)=0
ApproximationEnum
@ ApproximationEnum
Definition: EnumDefinitions.h:470
InversionCostFunctionsEnum
@ InversionCostFunctionsEnum
Definition: EnumDefinitions.h:211
Element::dViscositydBSSA
void dViscositydBSSA(IssmDouble *pdmudB, int dim, IssmDouble *xyz_list, Gauss *gauss, Input2 *vx_input, Input2 *vy_input)
Definition: Element.cpp:837
VzEnum
@ VzEnum
Definition: EnumDefinitions.h:853
EnumToStringx
const char * EnumToStringx(int enum_in)
Definition: EnumToStringx.cpp:15
SurfaceAreaEnum
@ SurfaceAreaEnum
Definition: EnumDefinitions.h:820
AdjointHorizAnalysis::GradientJBGradient
void GradientJBGradient(Element *element, Vector< IssmDouble > *gradient, int control_index)
Definition: AdjointHorizAnalysis.cpp:1357
DomainDimensionEnum
@ DomainDimensionEnum
Definition: EnumDefinitions.h:123
GsetEnum
@ GsetEnum
Definition: EnumDefinitions.h:1093
ThicknessAlongGradientEnum
@ ThicknessAlongGradientEnum
Definition: EnumDefinitions.h:839
alpha
IssmDouble alpha(IssmDouble x, IssmDouble y, IssmDouble z, int testid)
Definition: fsanalyticals.cpp:221
AdjointxEnum
@ AdjointxEnum
Definition: EnumDefinitions.h:466
Element::GetVerticesCoordinates
void GetVerticesCoordinates(IssmDouble **xyz_list)
Definition: Element.cpp:1446
AdjointHorizAnalysis::CreatePVectorL1L2
ElementVector * CreatePVectorL1L2(Element *element)
Definition: AdjointHorizAnalysis.cpp:584
AdjointHorizAnalysis::GradientJDragHydroSSA
void GradientJDragHydroSSA(Element *element, Vector< IssmDouble > *gradient, int control_index)
Definition: AdjointHorizAnalysis.cpp:2165
AdjointHorizAnalysis::GradientJBinitial
void GradientJBinitial(Element *element, Vector< IssmDouble > *gradient, int control_index)
Definition: AdjointHorizAnalysis.cpp:1423
AdjointHorizAnalysis::CreatePVectorSSA
ElementVector * CreatePVectorSSA(Element *element)
Definition: AdjointHorizAnalysis.cpp:824
AdjointHorizAnalysis::GradientJBbarSSA
void GradientJBbarSSA(Element *element, Vector< IssmDouble > *gradient, int control_index)
Definition: AdjointHorizAnalysis.cpp:1269
AdjointzEnum
@ AdjointzEnum
Definition: EnumDefinitions.h:468
Element::NumberofNodesPressure
virtual int NumberofNodesPressure(void)=0
AdjointHorizAnalysis::GradientJDSSA
void GradientJDSSA(Element *element, Vector< IssmDouble > *gradient, int control_index)
Definition: AdjointHorizAnalysis.cpp:2291
Domain3DEnum
@ Domain3DEnum
Definition: EnumDefinitions.h:536
AdjointHorizAnalysis::GradientJBFS
void GradientJBFS(Element *element, Vector< IssmDouble > *gradient, int control_index)
Definition: AdjointHorizAnalysis.cpp:1353
Friction
Definition: Friction.h:13
AdjointHorizAnalysis::CreatePVectorHO
ElementVector * CreatePVectorHO(Element *element)
Definition: AdjointHorizAnalysis.cpp:589
AdjointHorizAnalysis::CreatePVectorFS
ElementVector * CreatePVectorFS(Element *element)
Definition: AdjointHorizAnalysis.cpp:343
DamageDbarEnum
@ DamageDbarEnum
Definition: EnumDefinitions.h:518
Input2
Definition: Input2.h:18
Element::GetDofListLocal
void GetDofListLocal(int **pdoflist, int approximation_enum, int setenum)
Definition: Element.cpp:984
Element::NodalFunctionsVelocity
virtual void NodalFunctionsVelocity(IssmDouble *basis, Gauss *gauss)=0
LACrouzeixRaviartEnum
@ LACrouzeixRaviartEnum
Definition: EnumDefinitions.h:1138
Element::IsIceInElement
bool IsIceInElement()
Definition: Element.cpp:2021
Element::StrainRateHO
void StrainRateHO(IssmDouble *epsilon, IssmDouble *xyz_list, Gauss *gauss, Input2 *vx_input, Input2 *vy_input)
Definition: Element.cpp:4084
Element::GetVerticesCoordinatesBase
virtual void GetVerticesCoordinatesBase(IssmDouble **xyz_list)=0
_error_
#define _error_(StreamArgs)
Definition: exceptions.h:49
MaterialsRheologyBEnum
@ MaterialsRheologyBEnum
Definition: EnumDefinitions.h:643
RheologyBAbsGradientEnum
@ RheologyBAbsGradientEnum
Definition: EnumDefinitions.h:671
Gauss::begin
virtual int begin(void)=0
VxEnum
@ VxEnum
Definition: EnumDefinitions.h:846
AdjointHorizAnalysis::GradientJBbarGradient
void GradientJBbarGradient(Element *element, Vector< IssmDouble > *gradient, int control_index)
Definition: AdjointHorizAnalysis.cpp:1184
InversionNumCostFunctionsEnum
@ InversionNumCostFunctionsEnum
Definition: EnumDefinitions.h:224
SurfaceLogVelMisfitEnum
@ SurfaceLogVelMisfitEnum
Definition: EnumDefinitions.h:825
AdjointHorizAnalysis::CreateKMatrixFS
ElementMatrix * CreateKMatrixFS(Element *element)
Definition: AdjointHorizAnalysis.cpp:57
Element::JacobianDeterminant
virtual void JacobianDeterminant(IssmDouble *Jdet, IssmDouble *xyz_list, Gauss *gauss)=0
Gauss::GaussPoint
virtual void GaussPoint(int ig)=0
Element::GetInputValue
void GetInputValue(bool *pvalue, int enum_type)
Definition: Element.cpp:1177
Element::VelocityInterpolation
virtual int VelocityInterpolation()=0
ThicknessEnum
@ ThicknessEnum
Definition: EnumDefinitions.h:840
SSAApproximationEnum
@ SSAApproximationEnum
Definition: EnumDefinitions.h:1255
ElementVector
Definition: ElementVector.h:20
Element::StrainRateSSA
void StrainRateSSA(IssmDouble *epsilon, IssmDouble *xyz_list, Gauss *gauss, Input2 *vx_input, Input2 *vy_input)
Definition: Element.cpp:4138
Input2::GetInputValue
virtual void GetInputValue(IssmDouble *pvalue, Gauss *gauss)
Definition: Input2.h:38
InversionVxObsEnum
@ InversionVxObsEnum
Definition: EnumDefinitions.h:633
Gauss::weight
IssmDouble weight
Definition: Gauss.h:11
Element::GetVerticesCoordinatesTop
virtual void GetVerticesCoordinatesTop(IssmDouble **xyz_list)=0
SurfaceAverageVelMisfitEnum
@ SurfaceAverageVelMisfitEnum
Definition: EnumDefinitions.h:821
Element::IsOnSurface
bool IsOnSurface()
Definition: Element.cpp:1981
Element::NewGaussTop
virtual Gauss * NewGaussTop(int order)=0
Material::ViscositySSADerivativeEpsSquare
virtual void ViscositySSADerivativeEpsSquare(IssmDouble *pmu_prime, IssmDouble *epsilon, Gauss *gauss)=0
ThicknessAbsGradientEnum
@ ThicknessAbsGradientEnum
Definition: EnumDefinitions.h:836
Element::GetNumberOfVertices
virtual int GetNumberOfVertices(void)=0
AdjointHorizAnalysis::GradientJBSSA
void GradientJBSSA(Element *element, Vector< IssmDouble > *gradient, int control_index)
Definition: AdjointHorizAnalysis.cpp:1563
Domain2DverticalEnum
@ Domain2DverticalEnum
Definition: EnumDefinitions.h:535
ElementMatrix
Definition: ElementMatrix.h:19
FrictionCEnum
@ FrictionCEnum
Definition: EnumDefinitions.h:572
AdjointHorizAnalysis::GradientJDragHydroHO
void GradientJDragHydroHO(Element *element, Vector< IssmDouble > *gradient, int control_index)
Definition: AdjointHorizAnalysis.cpp:2092
Element::GradientIndexing
void GradientIndexing(int *indexing, int control_index)
Definition: Element.cpp:1510
RheologyBInitialguessMisfitEnum
@ RheologyBInitialguessMisfitEnum
Definition: EnumDefinitions.h:673
Gauss
Definition: Gauss.h:8
Element::material
Material * material
Definition: Element.h:50
Element::GetDofListLocalVelocity
void GetDofListLocalVelocity(int **pdoflist, int setenum)
Definition: Element.cpp:1078
Element::NodalFunctionsDerivatives
virtual void NodalFunctionsDerivatives(IssmDouble *dbasis, IssmDouble *xyz_list, Gauss *gauss)=0
RheologyBbarAbsGradientEnum
@ RheologyBbarAbsGradientEnum
Definition: EnumDefinitions.h:674
Element::NumberofNodesVelocity
virtual int NumberofNodesVelocity(void)=0
AdjointHorizAnalysis::GradientJDragSSA
void GradientJDragSSA(Element *element, Vector< IssmDouble > *gradient, int control_index)
Definition: AdjointHorizAnalysis.cpp:1904
AdjointHorizAnalysis::GradientJBbarHO
void GradientJBbarHO(Element *element, Vector< IssmDouble > *gradient, int control_index)
Definition: AdjointHorizAnalysis.cpp:1264
AdjointHorizAnalysis::GradientJBbarL1L2
void GradientJBbarL1L2(Element *element, Vector< IssmDouble > *gradient, int control_index)
Definition: AdjointHorizAnalysis.cpp:1259
ElementMatrix::values
IssmDouble * values
Definition: ElementMatrix.h:26
AdjointHorizAnalysis::CreateKMatrixSSA
ElementMatrix * CreateKMatrixSSA(Element *element)
Definition: AdjointHorizAnalysis.cpp:234
ThicknessAcrossGradientEnum
@ ThicknessAcrossGradientEnum
Definition: EnumDefinitions.h:838
Generated on Thu Jul 2 2020 08:09:26 for Ice Sheet System Model by   doxygen 1.8.19