Ice Sheet System Model  4.18
Code documentation
•All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros
Public Member Functions | Data Fields
Numericalflux Class Reference

#include <Numericalflux.h>

Inheritance diagram for Numericalflux:
Load Object

Public Member Functions

 Numericalflux ()
 
 Numericalflux (int numericalflux_id, int i, int index, IoModel *iomodel)
 
 ~Numericalflux ()
 
Object * copy ()
 
void DeepEcho ()
 
void Echo ()
 
int Id ()
 
void Marshall (char **pmarshalled_data, int *pmarshalled_data_size, int marshall_direction)
 
int ObjectEnum ()
 
void InputUpdateFromConstant (IssmDouble constant, int name)
 
void InputUpdateFromConstant (int constant, int name)
 
void InputUpdateFromConstant (bool constant, int name)
 
void InputUpdateFromIoModel (int index, IoModel *iomodel)
 
void InputUpdateFromMatrixDakota (IssmDouble *matrix, int nrows, int ncols, int name, int type)
 
void InputUpdateFromVector (IssmDouble *vector, int name, int type)
 
void InputUpdateFromVectorDakota (IssmDouble *vector, int name, int type)
 
void Configure (Elements *elements, Loads *loads, Nodes *nodes, Vertices *vertices, Materials *materials, Parameters *parameters)
 
void CreateJacobianMatrix (Matrix< IssmDouble > *Jff)
 
void CreateKMatrix (Matrix< IssmDouble > *Kff, Matrix< IssmDouble > *Kfs)
 
void CreatePVector (Vector< IssmDouble > *pf)
 
void GetNodesLidList (int *lidlist)
 
void GetNodesSidList (int *sidlist)
 
int GetNumberOfNodes (void)
 
int GetNumberOfNodesOneSide (void)
 
bool IsPenalty (void)
 
void PenaltyCreateJacobianMatrix (Matrix< IssmDouble > *Jff, IssmDouble kmax)
 
void PenaltyCreateKMatrix (Matrix< IssmDouble > *Kff, Matrix< IssmDouble > *kfs, IssmDouble kmax)
 
void PenaltyCreatePVector (Vector< IssmDouble > *pf, IssmDouble kmax)
 
void ResetHooks ()
 
void SetCurrentConfiguration (Elements *elements, Loads *loads, Nodes *nodes, Vertices *vertices, Materials *materials, Parameters *parameters)
 
void SetwiseNodeConnectivity (int *d_nz, int *o_nz, Node *node, bool *flags, int *flagsindices, int set1_enum, int set2_enum)
 
ElementMatrix * CreateKMatrixAdjointBalancethickness (void)
 
ElementMatrix * CreateKMatrixAdjointBalancethicknessBoundary (void)
 
ElementMatrix * CreateKMatrixAdjointBalancethicknessInternal (void)
 
ElementMatrix * CreateKMatrixBalancethickness (void)
 
ElementMatrix * CreateKMatrixBalancethicknessBoundary (void)
 
ElementMatrix * CreateKMatrixBalancethicknessInternal (void)
 
ElementMatrix * CreateKMatrixMasstransport (void)
 
ElementMatrix * CreateKMatrixMasstransportBoundary (void)
 
ElementMatrix * CreateKMatrixMasstransportInternal (void)
 
ElementVector * CreatePVectorAdjointBalancethickness (void)
 
ElementVector * CreatePVectorBalancethickness (void)
 
ElementVector * CreatePVectorBalancethicknessBoundary (void)
 
ElementVector * CreatePVectorBalancethicknessInternal (void)
 
ElementVector * CreatePVectorMasstransport (void)
 
ElementVector * CreatePVectorMasstransportBoundary (void)
 
ElementVector * CreatePVectorMasstransportInternal (void)
 
void GetNormal (IssmDouble *normal, IssmDouble xyz_list[4][3])
 
- Public Member Functions inherited from Load
virtual ~Load ()
 
- Public Member Functions inherited from Object
virtual ~Object ()
 

Data Fields

int id
 
int flux_type
 
int flux_degree
 
Hook * helement
 
Hook * hnodes
 
Hook * hvertices
 
Element * element
 
Vertex ** vertices
 
Node ** nodes
 
Parameters * parameters
 

Detailed Description

Definition at line 18 of file Numericalflux.h.

Constructor & Destructor Documentation

◆ Numericalflux() [1/2]

Numericalflux::Numericalflux ( )

Definition at line 21 of file Numericalflux.cpp.

21  {/*{{{*/
22  this->parameters = NULL;
23  this->helement = NULL;
24  this->element = NULL;
25  this->hnodes = NULL;
26  this->hvertices = NULL;
27  this->nodes = NULL;
28 }

◆ Numericalflux() [2/2]

Numericalflux::Numericalflux ( int  numericalflux_id,
int  i,
int  index,
IoModel *  iomodel 
)

Definition at line 30 of file Numericalflux.cpp.

30  {/*{{{*/
31 
32  /* Intermediary */
33  int pos1,pos2,pos3,pos4;
34  int numnodes;
35 
36  /*numericalflux constructor data: */
37  int numericalflux_elem_ids[2];
38  int numericalflux_vertex_ids[2];
39  int numericalflux_node_ids[4];
40  int numericalflux_type;
41  int numericalflux_degree;
42 
43  /*Get edge*/
44  int i1 = iomodel->faces[4*index+0];
45  int i2 = iomodel->faces[4*index+1];
46  int e1 = iomodel->faces[4*index+2];
47  int e2 = iomodel->faces[4*index+3];
48 
49  /*First, see wether this is an internal or boundary edge (if e2=-1)*/
50  if(e2==-1){
51  /* Boundary edge, only one element */
52  numericalflux_type=BoundaryEnum;
53  numericalflux_elem_ids[0]=e1;
54  }
55  else{
56  /* internal edge: connected to 2 elements */
57  numericalflux_type=InternalEnum;
58  numericalflux_elem_ids[0]=e1;
59  numericalflux_elem_ids[1]=e2;
60  }
61 
62  /*FIXME: hardcode element degree for now*/
63  this->flux_degree= P1DGEnum;
64  //this->flux_degree= P0DGEnum;
65 
66  /*1: Get vertices ids*/
67  numericalflux_vertex_ids[0]=i1;
68  numericalflux_vertex_ids[1]=i2;
69 
70  /*2: Get node ids*/
71  if(numericalflux_type==InternalEnum){
72  /*Get the column where these ids are located in the index*/
73  pos1=pos2=pos3=pos4=UNDEF;
74  for(int j=0;j<3;j++){
75  if(iomodel->elements[3*(e1-1)+j]==i1) pos1=j+1;
76  if(iomodel->elements[3*(e1-1)+j]==i2) pos2=j+1;
77  if(iomodel->elements[3*(e2-1)+j]==i1) pos3=j+1;
78  if(iomodel->elements[3*(e2-1)+j]==i2) pos4=j+1;
79  }
80  _assert_(pos1!=UNDEF && pos2!=UNDEF && pos3!=UNDEF && pos4!=UNDEF);
81 
82  /* We have the id of the elements and the position of the vertices in the index
83  * we can compute their dofs!*/
84  numericalflux_node_ids[0]=3*(e1-1)+pos1;
85  numericalflux_node_ids[1]=3*(e1-1)+pos2;
86  numericalflux_node_ids[2]=3*(e2-1)+pos3;
87  numericalflux_node_ids[3]=3*(e2-1)+pos4;
88  }
89  else{
90  /*Get the column where these ids are located in the index*/
91  pos1=pos2=UNDEF;
92  for(int j=0;j<3;j++){
93  if(iomodel->elements[3*(e1-1)+j]==i1) pos1=j+1;
94  if(iomodel->elements[3*(e1-1)+j]==i2) pos2=j+1;
95  }
96  _assert_(pos1!=UNDEF && pos2!=UNDEF);
97 
98  /* We have the id of the elements and the position of the vertices in the index
99  * we can compute their dofs!*/
100  numericalflux_node_ids[0]=3*(e1-1)+pos1;
101  numericalflux_node_ids[1]=3*(e1-1)+pos2;
102  }
103 
104  switch(this->flux_degree){
105  case P0DGEnum:
106  if(numericalflux_type==InternalEnum) numnodes = 2;
107  else numnodes = 1;
108  for(int i=0;i<numnodes;i++) numericalflux_node_ids[i] = numericalflux_elem_ids[i];
109  numericalflux_node_ids[1] = numericalflux_elem_ids[1];
110  break;
111  case P1DGEnum:
112  if(numericalflux_type==InternalEnum) numnodes = 4;
113  else numnodes = 2;
114  for(int i=0;i<numnodes;i++) numericalflux_node_ids[i] = numericalflux_node_ids[i]; //FIXME: to be improved...
115  break;
116  default:
117  _error_("not supported yet");
118 
119  }
120 
121  /*Assign object fields: */
122  this->id = numericalflux_id;
123  this->flux_type = numericalflux_type;
124 
125  /*Hooks: */
126  this->hnodes = new Hook(numericalflux_node_ids,numnodes);
127  this->hvertices = new Hook(&numericalflux_vertex_ids[0],2);
128  this->helement = new Hook(numericalflux_elem_ids,1); // take only the first element for now
129 
130  /*other fields*/
131  this->parameters = NULL;
132  this->element = NULL;
133  this->nodes = NULL;
134 }

◆ ~Numericalflux()

Numericalflux::~Numericalflux ( )

Definition at line 136 of file Numericalflux.cpp.

136  {/*{{{*/
137  this->parameters=NULL;
138  delete helement;
139  delete hnodes;
140  delete hvertices;
141 }

Member Function Documentation

◆ copy()

Object * Numericalflux::copy ( void  )
virtual

Implements Object.

Definition at line 145 of file Numericalflux.cpp.

145  {/*{{{*/
146 
147  Numericalflux* numericalflux=NULL;
148 
149  numericalflux=new Numericalflux();
150 
151  /*copy fields: */
152  numericalflux->id=this->id;
153  numericalflux->flux_type=this->flux_type;
154  numericalflux->flux_degree=this->flux_degree;
155 
156  /*point parameters: */
157  numericalflux->parameters=this->parameters;
158 
159  /*now deal with hooks and objects: */
160  numericalflux->hnodes = (Hook*)this->hnodes->copy();
161  numericalflux->hvertices = (Hook*)this->hvertices->copy();
162  numericalflux->helement = (Hook*)this->helement->copy();
163 
164  /*corresponding fields*/
165  numericalflux->nodes = (Node**)numericalflux->hnodes->deliverp();
166  numericalflux->vertices = (Vertex**)numericalflux->hvertices->deliverp();
167  numericalflux->element = (Element*)numericalflux->helement->delivers();
168 
169  return numericalflux;
170 }

◆ DeepEcho()

void Numericalflux::DeepEcho ( void  )
virtual

Implements Object.

Definition at line 172 of file Numericalflux.cpp.

172  {/*{{{*/
173 
174  _printf_("Numericalflux:\n");
175  _printf_(" id: " << id << "\n");
176  _printf_(" flux_type: " << this->flux_type<< "\n");
177  _printf_(" flux_degree: " << this->flux_degree<< "\n");
178  hnodes->DeepEcho();
179  hvertices->DeepEcho();
180  helement->DeepEcho();
181  _printf_(" parameters\n");
182  if(parameters)
183  parameters->DeepEcho();
184  else
185  _printf_(" NULL\n");
186 }

◆ Echo()

void Numericalflux::Echo ( void  )
virtual

Implements Object.

Definition at line 188 of file Numericalflux.cpp.

188  {/*{{{*/
189  _printf_("Numericalflux:\n");
190  _printf_(" id: " << id << "\n");
191  _printf_(" flux_type: " << this->flux_type<< "\n");
192  _printf_(" flux_degree: " << this->flux_degree<< "\n");
193  hnodes->Echo();
194  hvertices->Echo();
195  helement->Echo();
196  _printf_(" parameters: " << parameters << "\n");
197 }

◆ Id()

int Numericalflux::Id ( void  )
virtual

Implements Object.

Definition at line 199 of file Numericalflux.cpp.

199  {/*{{{*/
200  return id;
201 }

◆ Marshall()

void Numericalflux::Marshall ( char **  pmarshalled_data,
int *  pmarshalled_data_size,
int  marshall_direction 
)
virtual

Implements Object.

Definition at line 203 of file Numericalflux.cpp.

203  { /*{{{*/
204 
205  _assert_(this);
206 
207  /*ok, marshall operations: */
208  MARSHALLING_ENUM(NumericalfluxEnum);
209  MARSHALLING(id);
210  MARSHALLING(flux_type);
211  MARSHALLING(flux_degree);
212 
213  if(marshall_direction==MARSHALLING_BACKWARD){
214  this->hnodes = new Hook();
215  this->hvertices = new Hook();
216  this->helement = new Hook();
217  }
218 
219  this->hnodes->Marshall(pmarshalled_data,pmarshalled_data_size,marshall_direction);
220  this->helement->Marshall(pmarshalled_data,pmarshalled_data_size,marshall_direction);
221  this->hvertices->Marshall(pmarshalled_data,pmarshalled_data_size,marshall_direction);
222 
223  /*corresponding fields*/
224  nodes =(Node**)this->hnodes->deliverp();
225  vertices =(Vertex**)this->hvertices->deliverp();
226  element =(Element*)this->helement->delivers();
227 
228 }

◆ ObjectEnum()

int Numericalflux::ObjectEnum ( void  )
virtual

Implements Object.

Definition at line 230 of file Numericalflux.cpp.

230  {/*{{{*/
231 
232  return NumericalfluxEnum;
233 
234 }

◆ InputUpdateFromConstant() [1/3]

void Numericalflux::InputUpdateFromConstant ( IssmDouble  constant,
int  name 
)
inline

Definition at line 50 of file Numericalflux.h.

50 {/*Do nothing*/};

◆ InputUpdateFromConstant() [2/3]

void Numericalflux::InputUpdateFromConstant ( int  constant,
int  name 
)
inline

Definition at line 51 of file Numericalflux.h.

51 {/*Do nothing*/};

◆ InputUpdateFromConstant() [3/3]

void Numericalflux::InputUpdateFromConstant ( bool  constant,
int  name 
)
inline

Definition at line 52 of file Numericalflux.h.

52 {_error_("Not implemented yet!");}

◆ InputUpdateFromIoModel()

void Numericalflux::InputUpdateFromIoModel ( int  index,
IoModel *  iomodel 
)
inline

Definition at line 53 of file Numericalflux.h.

53 {_error_("not implemented yet");};

◆ InputUpdateFromMatrixDakota()

void Numericalflux::InputUpdateFromMatrixDakota ( IssmDouble *  matrix,
int  nrows,
int  ncols,
int  name,
int  type 
)
inline

Definition at line 54 of file Numericalflux.h.

54 {/*Do nothing*/}

◆ InputUpdateFromVector()

void Numericalflux::InputUpdateFromVector ( IssmDouble *  vector,
int  name,
int  type 
)
inline

Definition at line 55 of file Numericalflux.h.

55 {/*Do nothing*/}

◆ InputUpdateFromVectorDakota()

void Numericalflux::InputUpdateFromVectorDakota ( IssmDouble *  vector,
int  name,
int  type 
)
inline

Definition at line 56 of file Numericalflux.h.

56 {/*Do nothing*/}

◆ Configure()

void Numericalflux::Configure ( Elements *  elements,
Loads *  loads,
Nodes *  nodes,
Vertices *  vertices,
Materials *  materials,
Parameters *  parameters 
)
virtual

Implements Load.

Definition at line 238 of file Numericalflux.cpp.

238  {/*{{{*/
239 
240  /*Take care of hooking up all objects for this element, ie links the objects in the hooks to their respective
241  * datasets, using internal ids and offsets hidden in hooks: */
242  hnodes->configure((DataSet*)nodesin);
243  hvertices->configure((DataSet*)verticesin);
244  helement->configure((DataSet*)elementsin);
245 
246  /*Initialize hooked fields*/
247  this->nodes = (Node**)hnodes->deliverp();
248  this->vertices = (Vertex**)hvertices->deliverp();
249  this->element = (Element*)helement->delivers();
250 
251  /*point parameters to real dataset: */
252  this->parameters=parametersin;
253 }

◆ CreateJacobianMatrix()

void Numericalflux::CreateJacobianMatrix ( Matrix< IssmDouble > *  Jff)
inlinevirtual

Implements Load.

Definition at line 60 of file Numericalflux.h.

60 {_error_("Not implemented yet");};

◆ CreateKMatrix()

void Numericalflux::CreateKMatrix ( Matrix< IssmDouble > *  Kff,
Matrix< IssmDouble > *  Kfs 
)
virtual

Implements Load.

Definition at line 255 of file Numericalflux.cpp.

255  {/*{{{*/
256 
257  /*recover some parameters*/
258  ElementMatrix* Ke=NULL;
259  int analysis_type;
260  this->parameters->FindParam(&analysis_type,AnalysisTypeEnum);
261 
262  /*Just branch to the correct element stiffness matrix generator, according to the type of analysis we are carrying out: */
263  switch(analysis_type){
264  case MasstransportAnalysisEnum:
265  Ke=CreateKMatrixMasstransport();
266  break;
267  case BalancethicknessAnalysisEnum:
268  Ke=CreateKMatrixBalancethickness();
269  break;
270  case AdjointBalancethicknessAnalysisEnum:
271  Ke=CreateKMatrixAdjointBalancethickness();
272  break;
273  default:
274  _error_("analysis " << analysis_type << " (" << EnumToStringx(analysis_type) << ") not supported yet");
275  }
276 
277  /*Add to global matrix*/
278  if(Ke){
279  Ke->AddToGlobal(Kff,Kfs);
280  delete Ke;
281  }
282 
283 }

◆ CreatePVector()

void Numericalflux::CreatePVector ( Vector< IssmDouble > *  pf)
virtual

Implements Load.

Definition at line 285 of file Numericalflux.cpp.

285  {/*{{{*/
286 
287  /*recover some parameters*/
288  ElementVector* pe=NULL;
289  int analysis_type;
290  this->parameters->FindParam(&analysis_type,AnalysisTypeEnum);
291 
292  switch(analysis_type){
293  case MasstransportAnalysisEnum:
294  pe=CreatePVectorMasstransport();
295  break;
296  case BalancethicknessAnalysisEnum:
297  pe=CreatePVectorBalancethickness();
298  break;
299  case AdjointBalancethicknessAnalysisEnum:
300  pe=CreatePVectorAdjointBalancethickness();
301  break;
302  default:
303  _error_("analysis " << analysis_type << " (" << EnumToStringx(analysis_type) << ") not supported yet");
304  }
305 
306  /*Add to global matrix*/
307  if(pe){
308  pe->AddToGlobal(pf);
309  delete pe;
310  }
311 
312 }

◆ GetNodesLidList()

void Numericalflux::GetNodesLidList ( int *  lidlist)
virtual

Implements Load.

Definition at line 314 of file Numericalflux.cpp.

314  {/*{{{*/
315 
316  _assert_(lidlist);
317  _assert_(nodes);
318 
319  int numnodes = this->GetNumberOfNodes();
320  for(int i=0;i<numnodes;i++) lidlist[i]=nodes[i]->Lid();
321 }

◆ GetNodesSidList()

void Numericalflux::GetNodesSidList ( int *  sidlist)
virtual

Implements Load.

Definition at line 323 of file Numericalflux.cpp.

323  {/*{{{*/
324 
325  _assert_(sidlist);
326  _assert_(nodes);
327 
328  int numnodes = this->GetNumberOfNodes();
329  for(int i=0;i<numnodes;i++) sidlist[i]=nodes[i]->Sid();
330 }

◆ GetNumberOfNodes()

int Numericalflux::GetNumberOfNodes ( void  )
virtual

Implements Load.

Definition at line 332 of file Numericalflux.cpp.

332  {/*{{{*/
333 
334  if(this->flux_degree==P0DGEnum){
335  switch(this->flux_type){
336  case InternalEnum:
337  return 2;
338  case BoundaryEnum:
339  return 1;
340  default:
341  _error_("Numericalflux type " << EnumToStringx(this->flux_type) << " not supported yet");
342  }
343  }
344  else if(this->flux_degree==P1DGEnum){
345  switch(this->flux_type){
346  case InternalEnum:
347  return 4;
348  case BoundaryEnum:
349  return 2;
350  default:
351  _error_("Numericalflux type " << EnumToStringx(this->flux_type) << " not supported yet");
352  }
353  }
354  else{
355  _error_("Numericalflux " << EnumToStringx(this->flux_degree) << " not supported yet");
356  }
357 
358 }

◆ GetNumberOfNodesOneSide()

int Numericalflux::GetNumberOfNodesOneSide ( void  )

Definition at line 360 of file Numericalflux.cpp.

360  {/*{{{*/
361 
362  if(this->flux_degree==P0DGEnum){
363  return 1;
364  }
365  else if(this->flux_degree==P1DGEnum){
366  return 2;
367  }
368  else{
369  _error_("Numericalflux " << EnumToStringx(this->flux_degree) << " not supported yet");
370  }
371 
372 }

◆ IsPenalty()

bool Numericalflux::IsPenalty ( void  )
virtual

Implements Load.

Definition at line 374 of file Numericalflux.cpp.

374  {/*{{{*/
375  return false;
376 }

◆ PenaltyCreateJacobianMatrix()

void Numericalflux::PenaltyCreateJacobianMatrix ( Matrix< IssmDouble > *  Jff,
IssmDouble  kmax 
)
inlinevirtual

Implements Load.

Definition at line 68 of file Numericalflux.h.

68 {_error_("Not implemented yet");};

◆ PenaltyCreateKMatrix()

void Numericalflux::PenaltyCreateKMatrix ( Matrix< IssmDouble > *  Kff,
Matrix< IssmDouble > *  kfs,
IssmDouble  kmax 
)
virtual

Implements Load.

Definition at line 378 of file Numericalflux.cpp.

378  {/*{{{*/
379 
380  /*No stiffness loads applied, do nothing: */
381  return;
382 
383 }

◆ PenaltyCreatePVector()

void Numericalflux::PenaltyCreatePVector ( Vector< IssmDouble > *  pf,
IssmDouble  kmax 
)
virtual

Implements Load.

Definition at line 385 of file Numericalflux.cpp.

385  {/*{{{*/
386 
387  /*No penalty loads applied, do nothing: */
388  return;
389 
390 }

◆ ResetHooks()

void Numericalflux::ResetHooks ( )
virtual

Implements Load.

Definition at line 392 of file Numericalflux.cpp.

392  {/*{{{*/
393 
394  this->nodes = NULL;
395  this->vertices = NULL;
396  this->element = NULL;
397  this->parameters = NULL;
398 
399  /*Get Element type*/
400  this->hnodes->reset();
401  this->hvertices->reset();
402  this->helement->reset();
403 
404 }

◆ SetCurrentConfiguration()

void Numericalflux::SetCurrentConfiguration ( Elements *  elements,
Loads *  loads,
Nodes *  nodes,
Vertices *  vertices,
Materials *  materials,
Parameters *  parameters 
)
virtual

Implements Load.

Definition at line 406 of file Numericalflux.cpp.

406  {/*{{{*/
407  /*Nothing to do :)*/
408 
409 }

◆ SetwiseNodeConnectivity()

void Numericalflux::SetwiseNodeConnectivity ( int *  d_nz,
int *  o_nz,
Node *  node,
bool *  flags,
int *  flagsindices,
int  set1_enum,
int  set2_enum 
)
virtual

Implements Load.

Definition at line 411 of file Numericalflux.cpp.

411  {/*{{{*/
412 
413  /*Output */
414  int d_nz = 0;
415  int o_nz = 0;
416 
417  /*Loop over all nodes*/
418  for(int i=0;i<this->GetNumberOfNodes();i++){
419 
420  if(!flags[this->nodes[i]->Lid()]){
421 
422  /*flag current node so that no other element processes it*/
423  flags[this->nodes[i]->Lid()]=true;
424 
425  int counter=0;
426  while(flagsindices[counter]>=0) counter++;
427  flagsindices[counter]=this->nodes[i]->Lid();
428 
429  /*if node is clone, we have an off-diagonal non-zero, else it is a diagonal non-zero*/
430  switch(set2_enum){
431  case FsetEnum:
432  if(nodes[i]->fsize){
433  if(this->nodes[i]->IsClone())
434  o_nz += 1;
435  else
436  d_nz += 1;
437  }
438  break;
439  case GsetEnum:
440  if(nodes[i]->gsize){
441  if(this->nodes[i]->IsClone())
442  o_nz += 1;
443  else
444  d_nz += 1;
445  }
446  break;
447  case SsetEnum:
448  if(nodes[i]->ssize){
449  if(this->nodes[i]->IsClone())
450  o_nz += 1;
451  else
452  d_nz += 1;
453  }
454  break;
455  default: _error_("not supported");
456  }
457  }
458  }
459 
460  /*Assign output pointers: */
461  *pd_nz=d_nz;
462  *po_nz=o_nz;
463 }

◆ CreateKMatrixAdjointBalancethickness()

ElementMatrix * Numericalflux::CreateKMatrixAdjointBalancethickness ( void  )

Definition at line 467 of file Numericalflux.cpp.

467  {/*{{{*/
468 
469  switch(this->flux_type){
470  case InternalEnum:
471  return CreateKMatrixAdjointBalancethicknessInternal();
472  case BoundaryEnum:
473  return CreateKMatrixAdjointBalancethicknessBoundary();
474  default:
475  _error_("type not supported yet");
476  }
477 }

◆ CreateKMatrixAdjointBalancethicknessBoundary()

ElementMatrix * Numericalflux::CreateKMatrixAdjointBalancethicknessBoundary ( void  )

Definition at line 479 of file Numericalflux.cpp.

479  {/*{{{*/
480 
481  ElementMatrix* Ke=CreateKMatrixBalancethicknessBoundary();
482  if(Ke) Ke->Transpose();
483  return Ke;
484 }

◆ CreateKMatrixAdjointBalancethicknessInternal()

ElementMatrix * Numericalflux::CreateKMatrixAdjointBalancethicknessInternal ( void  )

Definition at line 486 of file Numericalflux.cpp.

486  {/*{{{*/
487 
488  ElementMatrix* Ke=CreateKMatrixBalancethicknessInternal();
489  if (Ke) Ke->Transpose();
490  return Ke;
491 }

◆ CreateKMatrixBalancethickness()

ElementMatrix * Numericalflux::CreateKMatrixBalancethickness ( void  )

Definition at line 493 of file Numericalflux.cpp.

493  {/*{{{*/
494 
495  switch(this->flux_type){
496  case InternalEnum:
497  return CreateKMatrixBalancethicknessInternal();
498  case BoundaryEnum:
499  return CreateKMatrixBalancethicknessBoundary();
500  default:
501  _error_("type not supported yet");
502  }
503 }

◆ CreateKMatrixBalancethicknessBoundary()

ElementMatrix * Numericalflux::CreateKMatrixBalancethicknessBoundary ( void  )

Definition at line 505 of file Numericalflux.cpp.

505  {/*{{{*/
506 
507  /*Initialize Element matrix and return if necessary*/
508  Tria* tria=(Tria*)element;
509  if(!tria->IsIceInElement()) return NULL;
510 
511  /* Intermediaries*/
512  IssmDouble DL,Jdet,vx,vy,mean_vx,mean_vy;
513  IssmDouble xyz_list[NUMVERTICES][3];
514  IssmDouble normal[2];
515 
516  /*Retrieve all inputs and parameters*/
517  GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
518  Input2* vxaverage_input=tria->GetInput2(VxEnum); _assert_(vxaverage_input);
519  Input2* vyaverage_input=tria->GetInput2(VyEnum); _assert_(vyaverage_input);
520  GetNormal(&normal[0],xyz_list);
521 
522  /*Check wether it is an inflow or outflow BC (0 is the middle of the segment)*/
523  int index1=tria->GetVertexIndex(vertices[0]);
524  int index2=tria->GetVertexIndex(vertices[1]);
525 
526  GaussTria* gauss=new GaussTria();
527  gauss->GaussEdgeCenter(index1,index2);
528  vxaverage_input->GetInputValue(&mean_vx,gauss);
529  vyaverage_input->GetInputValue(&mean_vy,gauss);
530  delete gauss;
531 
532  IssmDouble UdotN=mean_vx*normal[0]+mean_vy*normal[1];
533  if(UdotN<=0){
534  return NULL; /*(u,n)<0 -> inflow, PenaltyCreatePVector will take care of it*/
535  }
536 
537  /*Initialize Element vector and other vectors*/
538  int numnodes = this->GetNumberOfNodes();
539  ElementMatrix *Ke = new ElementMatrix(nodes,numnodes,this->parameters);
540  IssmDouble *basis = xNew<IssmDouble>(numnodes);
541 
542  /* Start looping on the number of gaussian points: */
543  gauss=new GaussTria(index1,index2,2);
544  for(int ig=gauss->begin();ig<gauss->end();ig++){
545 
546  gauss->GaussPoint(ig);
547 
548  tria->GetSegmentNodalFunctions(&basis[0],gauss,index1,index2,tria->FiniteElement());
549 
550  vxaverage_input->GetInputValue(&vx,gauss);
551  vyaverage_input->GetInputValue(&vy,gauss);
552  UdotN=vx*normal[0]+vy*normal[1];
553  tria->GetSegmentJacobianDeterminant(&Jdet,&xyz_list[0][0],gauss);
554  DL=gauss->weight*Jdet*UdotN;
555 
556  for(int i=0;i<numnodes;i++){
557  for(int j=0;j<numnodes;j++){
558  Ke->values[i*numnodes+j]+=DL*basis[i]*basis[j];
559  }
560  }
561  }
562 
563  /*Clean up and return*/
564  xDelete<IssmDouble>(basis);
565  delete gauss;
566  return Ke;
567 }

◆ CreateKMatrixBalancethicknessInternal()

ElementMatrix * Numericalflux::CreateKMatrixBalancethicknessInternal ( void  )

Definition at line 569 of file Numericalflux.cpp.

569  {/*{{{*/
570 
571  /*Initialize Element matrix and return if necessary*/
572  Tria* tria=(Tria*)element;
573  if(!tria->IsIceInElement()) return NULL;
574 
575  /* Intermediaries*/
576  IssmDouble A1,A2,Jdet,vx,vy,UdotN;
577  IssmDouble xyz_list[NUMVERTICES][3];
578  IssmDouble normal[2];
579 
580  /*Fetch number of nodes for this flux*/
581  int numnodes = this->GetNumberOfNodes();
582  int numnodes_plus = this->GetNumberOfNodesOneSide();
583  int numnodes_minus = numnodes_plus; /*For now we are not doing p-adaptive DG*/
584  _assert_(numnodes==numnodes_plus+numnodes_minus);
585 
586  /*Initialize variables*/
587  ElementMatrix *Ke = new ElementMatrix(nodes,numnodes,this->parameters);
588  IssmDouble *basis_plus = xNew<IssmDouble>(numnodes_plus);
589  IssmDouble *basis_minus = xNew<IssmDouble>(numnodes_minus);
590 
591  /*Retrieve all inputs and parameters*/
592  GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
593  Input2* vxaverage_input=tria->GetInput2(VxEnum); _assert_(vxaverage_input);
594  Input2* vyaverage_input=tria->GetInput2(VyEnum); _assert_(vyaverage_input);
595  GetNormal(&normal[0],xyz_list);
596 
597  /* Start looping on the number of gaussian points: */
598  int index1=tria->GetVertexIndex(vertices[0]);
599  int index2=tria->GetVertexIndex(vertices[1]);
600  GaussTria* gauss=new GaussTria(index1,index2,2);
601  for(int ig=gauss->begin();ig<gauss->end();ig++){
602 
603  gauss->GaussPoint(ig);
604 
605  tria->GetSegmentNodalFunctions(&basis_plus[0] ,gauss,index1,index2,tria->FiniteElement());
606  tria->GetSegmentNodalFunctions(&basis_minus[0],gauss,index1,index2,tria->FiniteElement());
607 
608  vxaverage_input->GetInputValue(&vx,gauss);
609  vyaverage_input->GetInputValue(&vy,gauss);
610  UdotN=vx*normal[0]+vy*normal[1];
611  tria->GetSegmentJacobianDeterminant(&Jdet,&xyz_list[0][0],gauss);
612  A1=gauss->weight*Jdet*UdotN/2;
613  A2=gauss->weight*Jdet*fabs(UdotN)/2;
614 
615  /*Term 1 (numerical flux): {Hv}.[[phi]] = 0.5(H+v+ + H-v-)(phi+n+ + phi-n-)
616  * = v.n/2 (H+phi+ + H-phi+ -H+phi- -H-phi-)
617  * = v.n/2 (H+phi+ + H-phi+ -H+phi- -H-phi-)
618  *
619  *Term 2 (stabilization) |v.n|/2 [[H]].[[phi]] = |v.n|/2 (H+n+ + H-n-)(phi+n+ + phi-n-)
620  * = |v.n|/2 (H+phi+ -H-phi+ -H+phi- +H-phi-)
621  * | A++ | A+- |
622  * K = |-----------|
623  * | A-+ | A-- |
624  *
625  *These 4 terms for each expressions are added independently*/
626 
627  /*First term A++*/
628  for(int i=0;i<numnodes_plus;i++){
629  for(int j=0;j<numnodes_plus;j++){
630  Ke->values[i*numnodes+j] += A1*(basis_plus[j]*basis_plus[i]);
631  Ke->values[i*numnodes+j] += A2*(basis_plus[j]*basis_plus[i]);
632  }
633  }
634  /*Second term A+-*/
635  for(int i=0;i<numnodes_plus;i++){
636  for(int j=0;j<numnodes_minus;j++){
637  Ke->values[i*numnodes+numnodes_plus+j] += A1*(basis_minus[j]*basis_plus[i]);
638  Ke->values[i*numnodes+numnodes_plus+j] += -A2*(basis_minus[j]*basis_plus[i]);
639  }
640  }
641  /*Third term A-+*/
642  for(int i=0;i<numnodes_minus;i++){
643  for(int j=0;j<numnodes_plus;j++){
644  Ke->values[(numnodes_plus+i)*numnodes+j] += -A1*(basis_plus[j]*basis_minus[i]);
645  Ke->values[(numnodes_plus+i)*numnodes+j] += -A2*(basis_plus[j]*basis_minus[i]);
646  }
647  }
648  /*Fourth term A-+*/
649  for(int i=0;i<numnodes_minus;i++){
650  for(int j=0;j<numnodes_minus;j++){
651  Ke->values[(numnodes_plus+i)*numnodes+numnodes_plus+j] += -A1*(basis_minus[j]*basis_minus[i]);
652  Ke->values[(numnodes_plus+i)*numnodes+numnodes_plus+j] += A2*(basis_minus[j]*basis_minus[i]);
653  }
654  }
655  }
656 
657  /*Clean up and return*/
658  xDelete<IssmDouble>(basis_plus);
659  xDelete<IssmDouble>(basis_minus);
660  delete gauss;
661  return Ke;
662 }

◆ CreateKMatrixMasstransport()

ElementMatrix * Numericalflux::CreateKMatrixMasstransport ( void  )

Definition at line 664 of file Numericalflux.cpp.

664  {/*{{{*/
665 
666  switch(this->flux_type){
667  case InternalEnum:
668  return CreateKMatrixMasstransportInternal();
669  case BoundaryEnum:
670  return CreateKMatrixMasstransportBoundary();
671  default:
672  _error_("type not supported yet");
673  }
674 }

◆ CreateKMatrixMasstransportBoundary()

ElementMatrix * Numericalflux::CreateKMatrixMasstransportBoundary ( void  )

Definition at line 676 of file Numericalflux.cpp.

676  {/*{{{*/
677 
678  /*Initialize Element matrix and return if necessary*/
679  Tria* tria=(Tria*)element;
680  if(!tria->IsIceInElement()) return NULL;
681 
682  /* Intermediaries*/
683  IssmDouble DL,Jdet,vx,vy,mean_vx,mean_vy;
684  IssmDouble xyz_list[NUMVERTICES][3];
685  IssmDouble normal[2];
686 
687  /*Retrieve all inputs and parameters*/
688  GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
689  IssmDouble dt = parameters->FindParam(TimesteppingTimeStepEnum);
690  Input2* vxaverage_input=tria->GetInput2(VxEnum); _assert_(vxaverage_input);
691  Input2* vyaverage_input=tria->GetInput2(VyEnum); _assert_(vyaverage_input);
692  GetNormal(&normal[0],xyz_list);
693 
694  /*Check wether it is an inflow or outflow BC (0 is the middle of the segment)*/
695  int index1=tria->GetVertexIndex(vertices[0]);
696  int index2=tria->GetVertexIndex(vertices[1]);
697 
698  GaussTria* gauss=new GaussTria();
699  gauss->GaussEdgeCenter(index1,index2);
700  vxaverage_input->GetInputValue(&mean_vx,gauss);
701  vyaverage_input->GetInputValue(&mean_vy,gauss);
702  delete gauss;
703 
704  IssmDouble UdotN=mean_vx*normal[0]+mean_vy*normal[1];
705  if(UdotN<=0){
706  return NULL; /*(u,n)<0 -> inflow, PenaltyCreatePVector will take care of it*/
707  }
708 
709  /*Initialize Element vector and other vectors*/
710  int numnodes = this->GetNumberOfNodes();
711  ElementMatrix *Ke = new ElementMatrix(nodes,numnodes,this->parameters);
712  IssmDouble *basis = xNew<IssmDouble>(numnodes);
713 
714  /* Start looping on the number of gaussian points: */
715  gauss=new GaussTria(index1,index2,2);
716  for(int ig=gauss->begin();ig<gauss->end();ig++){
717 
718  gauss->GaussPoint(ig);
719 
720  tria->GetSegmentNodalFunctions(&basis[0],gauss,index1,index2,tria->FiniteElement());
721 
722  vxaverage_input->GetInputValue(&vx,gauss);
723  vyaverage_input->GetInputValue(&vy,gauss);
724  UdotN=vx*normal[0]+vy*normal[1];
725  tria->GetSegmentJacobianDeterminant(&Jdet,&xyz_list[0][0],gauss);
726  DL=gauss->weight*Jdet*dt*UdotN;
727 
728  for(int i=0;i<numnodes;i++){
729  for(int j=0;j<numnodes;j++){
730  Ke->values[i*numnodes+j]+=DL*basis[i]*basis[j];
731  }
732  }
733  }
734 
735  /*Clean up and return*/
736  xDelete<IssmDouble>(basis);
737  delete gauss;
738  return Ke;
739 }

◆ CreateKMatrixMasstransportInternal()

ElementMatrix * Numericalflux::CreateKMatrixMasstransportInternal ( void  )

Definition at line 741 of file Numericalflux.cpp.

741  {/*{{{*/
742 
743  /*Initialize Element matrix and return if necessary*/
744  Tria* tria=(Tria*)element;
745  if(!tria->IsIceInElement()) return NULL;
746 
747  /* Intermediaries*/
748  IssmDouble A1,A2,Jdet,vx,vy,UdotN;
749  IssmDouble xyz_list[NUMVERTICES][3];
750  IssmDouble normal[2];
751 
752  /*Fetch number of nodes for this flux*/
753  int numnodes = this->GetNumberOfNodes();
754  int numnodes_plus = this->GetNumberOfNodesOneSide();
755  int numnodes_minus = numnodes_plus; /*For now we are not doing p-adaptive DG*/
756  _assert_(numnodes==numnodes_plus+numnodes_minus);
757 
758  /*Initialize variables*/
759  ElementMatrix *Ke = new ElementMatrix(nodes,numnodes,this->parameters);
760  IssmDouble *basis_plus = xNew<IssmDouble>(numnodes_plus);
761  IssmDouble *basis_minus = xNew<IssmDouble>(numnodes_minus);
762 
763  /*Retrieve all inputs and parameters*/
764  GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
765  IssmDouble dt = parameters->FindParam(TimesteppingTimeStepEnum);
766  Input2* vxaverage_input=tria->GetInput2(VxEnum); _assert_(vxaverage_input);
767  Input2* vyaverage_input=tria->GetInput2(VyEnum); _assert_(vyaverage_input);
768  GetNormal(&normal[0],xyz_list);
769 
770  /* Start looping on the number of gaussian points: */
771  int index1=tria->GetVertexIndex(vertices[0]);
772  int index2=tria->GetVertexIndex(vertices[1]);
773  GaussTria* gauss=new GaussTria(index1,index2,2);
774  for(int ig=gauss->begin();ig<gauss->end();ig++){
775 
776  gauss->GaussPoint(ig);
777 
778  tria->GetSegmentNodalFunctions(&basis_plus[0] ,gauss,index1,index2,tria->FiniteElement());
779  tria->GetSegmentNodalFunctions(&basis_minus[0],gauss,index1,index2,tria->FiniteElement());
780 
781  vxaverage_input->GetInputValue(&vx,gauss);
782  vyaverage_input->GetInputValue(&vy,gauss);
783  UdotN=vx*normal[0]+vy*normal[1];
784  tria->GetSegmentJacobianDeterminant(&Jdet,&xyz_list[0][0],gauss);
785  A1=gauss->weight*Jdet*dt*UdotN/2;
786  A2=gauss->weight*Jdet*dt*fabs(UdotN)/2;
787 
788  /*Term 1 (numerical flux): {Hv}.[[phi]] = 0.5(H+v+ + H-v-)(phi+n+ + phi-n-)
789  * = v.n/2 (H+phi+ + H-phi+ -H+phi- -H-phi-)
790  * = v.n/2 (H+phi+ + H-phi+ -H+phi- -H-phi-)
791  *
792  *Term 2 (stabilization) |v.n|/2 [[H]].[[phi]] = |v.n|/2 (H+n+ + H-n-)(phi+n+ + phi-n-)
793  * = |v.n|/2 (H+phi+ -H-phi+ -H+phi- +H-phi-)
794  * | A++ | A+- |
795  * K = |-----------|
796  * | A-+ | A-- |
797  *
798  *These 4 terms for each expressions are added independently*/
799 
800  /*First term A++*/
801  for(int i=0;i<numnodes_plus;i++){
802  for(int j=0;j<numnodes_plus;j++){
803  Ke->values[i*numnodes+j] += A1*(basis_plus[j]*basis_plus[i]);
804  Ke->values[i*numnodes+j] += A2*(basis_plus[j]*basis_plus[i]);
805  }
806  }
807  /*Second term A+-*/
808  for(int i=0;i<numnodes_plus;i++){
809  for(int j=0;j<numnodes_minus;j++){
810  Ke->values[i*numnodes+numnodes_plus+j] += A1*(basis_minus[j]*basis_plus[i]);
811  Ke->values[i*numnodes+numnodes_plus+j] += -A2*(basis_minus[j]*basis_plus[i]);
812  }
813  }
814  /*Third term A-+*/
815  for(int i=0;i<numnodes_minus;i++){
816  for(int j=0;j<numnodes_plus;j++){
817  Ke->values[(numnodes_plus+i)*numnodes+j] += -A1*(basis_plus[j]*basis_minus[i]);
818  Ke->values[(numnodes_plus+i)*numnodes+j] += -A2*(basis_plus[j]*basis_minus[i]);
819  }
820  }
821  /*Fourth term A-+*/
822  for(int i=0;i<numnodes_minus;i++){
823  for(int j=0;j<numnodes_minus;j++){
824  Ke->values[(numnodes_plus+i)*numnodes+numnodes_plus+j] += -A1*(basis_minus[j]*basis_minus[i]);
825  Ke->values[(numnodes_plus+i)*numnodes+numnodes_plus+j] += A2*(basis_minus[j]*basis_minus[i]);
826  }
827  }
828  }
829 
830  /*Clean up and return*/
831  xDelete<IssmDouble>(basis_plus);
832  xDelete<IssmDouble>(basis_minus);
833  delete gauss;
834  return Ke;
835 }

◆ CreatePVectorAdjointBalancethickness()

ElementVector * Numericalflux::CreatePVectorAdjointBalancethickness ( void  )

Definition at line 837 of file Numericalflux.cpp.

837  {/*{{{*/
838 
839  /*No PVector for the Adjoint*/
840  return NULL;
841 }

◆ CreatePVectorBalancethickness()

ElementVector * Numericalflux::CreatePVectorBalancethickness ( void  )

Definition at line 843 of file Numericalflux.cpp.

843  {/*{{{*/
844 
845  switch(this->flux_type){
846  case InternalEnum:
847  return CreatePVectorBalancethicknessInternal();
848  case BoundaryEnum:
849  return CreatePVectorBalancethicknessBoundary();
850  default:
851  _error_("type not supported yet");
852  }
853 }

◆ CreatePVectorBalancethicknessBoundary()

ElementVector * Numericalflux::CreatePVectorBalancethicknessBoundary ( void  )

Definition at line 855 of file Numericalflux.cpp.

855  {/*{{{*/
856 
857  /*Initialize Load Vector and return if necessary*/
858  Tria* tria=(Tria*)element;
859  if(!tria->IsIceInElement()) return NULL;
860 
861  /* Intermediaries*/
862  IssmDouble DL,Jdet,vx,vy,mean_vx,mean_vy,thickness;
863  IssmDouble xyz_list[NUMVERTICES][3];
864  IssmDouble normal[2];
865 
866  /*Retrieve all inputs and parameters*/
867  GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
868  Input2* vxaverage_input = tria->GetInput2(VxEnum); _assert_(vxaverage_input);
869  Input2* vyaverage_input = tria->GetInput2(VyEnum); _assert_(vyaverage_input);
870  Input2* thickness_input = tria->GetInput2(ThicknessEnum); _assert_(thickness_input);
871  GetNormal(&normal[0],xyz_list);
872 
873  /*Check wether it is an inflow or outflow BC (0 is the middle of the segment)*/
874  int index1=tria->GetVertexIndex(vertices[0]);
875  int index2=tria->GetVertexIndex(vertices[1]);
876  GaussTria* gauss=new GaussTria();
877  gauss->GaussEdgeCenter(index1,index2);
878  vxaverage_input->GetInputValue(&mean_vx,gauss);
879  vyaverage_input->GetInputValue(&mean_vy,gauss);
880  delete gauss;
881  IssmDouble UdotN=mean_vx*normal[0]+mean_vy*normal[1];
882  if(UdotN>0){
883  return NULL; /*(u,n)>0 -> outflow, PenaltyCreateKMatrix will take care of it*/
884  }
885 
886  /*Initialize Load Vector */
887  int numnodes = this->GetNumberOfNodes();
888  ElementVector *pe = new ElementVector(nodes,numnodes,this->parameters);
889  IssmDouble *basis = xNew<IssmDouble>(numnodes);
890 
891  /* Start looping on the number of gaussian points: */
892  gauss=new GaussTria(index1,index2,2);
893  for(int ig=gauss->begin();ig<gauss->end();ig++){
894 
895  gauss->GaussPoint(ig);
896 
897  tria->GetSegmentNodalFunctions(&basis[0],gauss,index1,index2,tria->FiniteElement());
898 
899  vxaverage_input->GetInputValue(&vx,gauss);
900  vyaverage_input->GetInputValue(&vy,gauss);
901  thickness_input->GetInputValue(&thickness,gauss);
902 
903  UdotN=vx*normal[0]+vy*normal[1];
904  tria->GetSegmentJacobianDeterminant(&Jdet,&xyz_list[0][0],gauss);
905  DL= - gauss->weight*Jdet*UdotN*thickness;
906 
907  for(int i=0;i<numnodes;i++) pe->values[i] += DL*basis[i];
908  }
909 
910  /*Clean up and return*/
911  xDelete<IssmDouble>(basis);
912  delete gauss;
913  return pe;
914 }

◆ CreatePVectorBalancethicknessInternal()

ElementVector * Numericalflux::CreatePVectorBalancethicknessInternal ( void  )

Definition at line 916 of file Numericalflux.cpp.

916  {/*{{{*/
917 
918  /*Nothing added to PVector*/
919  return NULL;
920 
921 }

◆ CreatePVectorMasstransport()

ElementVector * Numericalflux::CreatePVectorMasstransport ( void  )

Definition at line 923 of file Numericalflux.cpp.

923  {/*{{{*/
924 
925  switch(this->flux_type){
926  case InternalEnum:
927  return CreatePVectorMasstransportInternal();
928  case BoundaryEnum:
929  return CreatePVectorMasstransportBoundary();
930  default:
931  _error_("type not supported yet");
932  }
933 }

◆ CreatePVectorMasstransportBoundary()

ElementVector * Numericalflux::CreatePVectorMasstransportBoundary ( void  )

Definition at line 935 of file Numericalflux.cpp.

935  {/*{{{*/
936 
937  /*Initialize Load Vector and return if necessary*/
938  Tria* tria=(Tria*)element;
939  if(!tria->IsIceInElement()) return NULL;
940 
941  /* Intermediaries*/
942  IssmDouble DL,Jdet,vx,vy,mean_vx,mean_vy,thickness;
943  IssmDouble xyz_list[NUMVERTICES][3];
944  IssmDouble normal[2];
945 
946  /*Retrieve all inputs and parameters*/
947  GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
948  IssmDouble dt = parameters->FindParam(TimesteppingTimeStepEnum);
949  Input2* vxaverage_input = tria->GetInput2(VxEnum); _assert_(vxaverage_input);
950  Input2* vyaverage_input = tria->GetInput2(VyEnum); _assert_(vyaverage_input);
951  Input2* spcthickness_input = tria->GetInput2(MasstransportSpcthicknessEnum); _assert_(spcthickness_input);
952  GetNormal(&normal[0],xyz_list);
953 
954  /*Check wether it is an inflow or outflow BC (0 is the middle of the segment)*/
955  int index1=tria->GetVertexIndex(vertices[0]);
956  int index2=tria->GetVertexIndex(vertices[1]);
957  GaussTria* gauss=new GaussTria();
958  gauss->GaussEdgeCenter(index1,index2);
959  vxaverage_input->GetInputValue(&mean_vx,gauss);
960  vyaverage_input->GetInputValue(&mean_vy,gauss);
961  delete gauss;
962  IssmDouble UdotN=mean_vx*normal[0]+mean_vy*normal[1];
963  if(UdotN>0){
964  return NULL; /*(u,n)>0 -> outflow, PenaltyCreateKMatrix will take care of it*/
965  }
966 
967  /*Initialize Load Vector */
968  int numnodes = this->GetNumberOfNodes();
969  ElementVector *pe = new ElementVector(nodes,numnodes,this->parameters);
970  IssmDouble *basis = xNew<IssmDouble>(numnodes);
971 
972  /* Start looping on the number of gaussian points: */
973  gauss=new GaussTria(index1,index2,2);
974  for(int ig=gauss->begin();ig<gauss->end();ig++){
975 
976  gauss->GaussPoint(ig);
977 
978  tria->GetSegmentNodalFunctions(&basis[0],gauss,index1,index2,tria->FiniteElement());
979 
980  vxaverage_input->GetInputValue(&vx,gauss);
981  vyaverage_input->GetInputValue(&vy,gauss);
982  spcthickness_input->GetInputValue(&thickness,gauss);
983  if(xIsNan<IssmDouble>(thickness)) _error_("Cannot weakly apply constraint because NaN was provided");
984 
985  UdotN=vx*normal[0]+vy*normal[1];
986  tria->GetSegmentJacobianDeterminant(&Jdet,&xyz_list[0][0],gauss);
987  DL= - gauss->weight*Jdet*dt*UdotN*thickness;
988 
989  for(int i=0;i<numnodes;i++) pe->values[i] += DL*basis[i];
990  }
991 
992  /*Clean up and return*/
993  xDelete<IssmDouble>(basis);
994  delete gauss;
995  return pe;
996 }

◆ CreatePVectorMasstransportInternal()

ElementVector * Numericalflux::CreatePVectorMasstransportInternal ( void  )

Definition at line 998 of file Numericalflux.cpp.

998  {/*{{{*/
999 
1000  /*Nothing added to PVector*/
1001  return NULL;
1002 
1003 }

◆ GetNormal()

void Numericalflux::GetNormal ( IssmDouble *  normal,
IssmDouble  xyz_list[4][3] 
)

Definition at line 1005 of file Numericalflux.cpp.

1005  {/*{{{*/
1006 
1007  /*Build unit outward pointing vector*/
1008  IssmDouble vector[2];
1009 
1010  vector[0]=xyz_list[1][0] - xyz_list[0][0];
1011  vector[1]=xyz_list[1][1] - xyz_list[0][1];
1012 
1013  IssmDouble norm=sqrt(pow(vector[0],2.0)+pow(vector[1],2.0));
1014 
1015  normal[0]= + vector[1]/norm;
1016  normal[1]= - vector[0]/norm;
1017 }

Field Documentation

◆ id

int Numericalflux::id

Definition at line 21 of file Numericalflux.h.

◆ flux_type

int Numericalflux::flux_type

Definition at line 22 of file Numericalflux.h.

◆ flux_degree

int Numericalflux::flux_degree

Definition at line 23 of file Numericalflux.h.

◆ helement

Hook* Numericalflux::helement

Definition at line 26 of file Numericalflux.h.

◆ hnodes

Hook* Numericalflux::hnodes

Definition at line 27 of file Numericalflux.h.

◆ hvertices

Hook* Numericalflux::hvertices

Definition at line 28 of file Numericalflux.h.

◆ element

Element* Numericalflux::element

Definition at line 31 of file Numericalflux.h.

◆ vertices

Vertex** Numericalflux::vertices

Definition at line 32 of file Numericalflux.h.

◆ nodes

Node** Numericalflux::nodes

Definition at line 33 of file Numericalflux.h.

◆ parameters

Parameters* Numericalflux::parameters

Definition at line 34 of file Numericalflux.h.

The documentation for this class was generated from the following files:
BalancethicknessAnalysisEnum
@ BalancethicknessAnalysisEnum
Definition: EnumDefinitions.h:981
Numericalflux::helement
Hook * helement
Definition: Numericalflux.h:26
_assert_
#define _assert_(ignore)
Definition: exceptions.h:37
IssmDouble
double IssmDouble
Definition: types.h:37
Numericalflux::CreateKMatrixBalancethicknessInternal
ElementMatrix * CreateKMatrixBalancethicknessInternal(void)
Definition: Numericalflux.cpp:569
MasstransportAnalysisEnum
@ MasstransportAnalysisEnum
Definition: EnumDefinitions.h:1163
Numericalflux::flux_degree
int flux_degree
Definition: Numericalflux.h:23
TriaRef::GetSegmentJacobianDeterminant
void GetSegmentJacobianDeterminant(IssmDouble *Jdet, IssmDouble *xyz_list, Gauss *gauss)
Definition: TriaRef.cpp:229
Numericalflux::CreatePVectorMasstransportBoundary
ElementVector * CreatePVectorMasstransportBoundary(void)
Definition: Numericalflux.cpp:935
Numericalflux::CreatePVectorAdjointBalancethickness
ElementVector * CreatePVectorAdjointBalancethickness(void)
Definition: Numericalflux.cpp:837
_printf_
#define _printf_(StreamArgs)
Definition: Print.h:22
Numericalflux::hnodes
Hook * hnodes
Definition: Numericalflux.h:27
Hook::deliverp
Object ** deliverp(void)
Definition: Hook.cpp:187
MARSHALLING_ENUM
#define MARSHALLING_ENUM(EN)
Definition: Marshalling.h:14
TimesteppingTimeStepEnum
@ TimesteppingTimeStepEnum
Definition: EnumDefinitions.h:433
Hook::DeepEcho
void DeepEcho(void)
Definition: Hook.cpp:77
Numericalflux::GetNormal
void GetNormal(IssmDouble *normal, IssmDouble xyz_list[4][3])
Definition: Numericalflux.cpp:1005
Numericalflux::GetNumberOfNodes
int GetNumberOfNodes(void)
Definition: Numericalflux.cpp:332
Numericalflux::flux_type
int flux_type
Definition: Numericalflux.h:22
SsetEnum
@ SsetEnum
Definition: EnumDefinitions.h:1282
BoundaryEnum
@ BoundaryEnum
Definition: EnumDefinitions.h:999
P1DGEnum
@ P1DGEnum
Definition: EnumDefinitions.h:1215
ElementVector::values
IssmDouble * values
Definition: ElementVector.h:24
GaussTria::begin
int begin(void)
Definition: GaussTria.cpp:356
Numericalflux::CreatePVectorMasstransportInternal
ElementVector * CreatePVectorMasstransportInternal(void)
Definition: Numericalflux.cpp:998
Numericalflux::element
Element * element
Definition: Numericalflux.h:31
Numericalflux::CreateKMatrixMasstransportInternal
ElementMatrix * CreateKMatrixMasstransportInternal(void)
Definition: Numericalflux.cpp:741
VyEnum
@ VyEnum
Definition: EnumDefinitions.h:850
TriaRef::GetSegmentNodalFunctions
void GetSegmentNodalFunctions(IssmDouble *basis, Gauss *gauss, int index1, int index2, int finiteelement)
Definition: TriaRef.cpp:243
Hook::reset
void reset(void)
Definition: Hook.cpp:211
Numericalflux::Numericalflux
Numericalflux()
Definition: Numericalflux.cpp:21
P0DGEnum
@ P0DGEnum
Definition: EnumDefinitions.h:1214
Element
Definition: Element.h:41
GaussTria
Definition: GaussTria.h:12
AdjointBalancethicknessAnalysisEnum
@ AdjointBalancethicknessAnalysisEnum
Definition: EnumDefinitions.h:971
Parameters::DeepEcho
void DeepEcho()
Definition: Parameters.cpp:99
MasstransportSpcthicknessEnum
@ MasstransportSpcthicknessEnum
Definition: EnumDefinitions.h:642
Hook::delivers
Object * delivers(void)
Definition: Hook.cpp:191
Numericalflux::CreateKMatrixAdjointBalancethicknessInternal
ElementMatrix * CreateKMatrixAdjointBalancethicknessInternal(void)
Definition: Numericalflux.cpp:486
Numericalflux::CreateKMatrixAdjointBalancethicknessBoundary
ElementMatrix * CreateKMatrixAdjointBalancethicknessBoundary(void)
Definition: Numericalflux.cpp:479
Numericalflux::id
int id
Definition: Numericalflux.h:21
NUMVERTICES
#define NUMVERTICES
Definition: Numericalflux.cpp:18
NumericalfluxEnum
@ NumericalfluxEnum
Definition: EnumDefinitions.h:1206
Numericalflux::CreateKMatrixBalancethickness
ElementMatrix * CreateKMatrixBalancethickness(void)
Definition: Numericalflux.cpp:493
EnumToStringx
const char * EnumToStringx(int enum_in)
Definition: EnumToStringx.cpp:15
Hook
Definition: Hook.h:16
Numericalflux::CreateKMatrixMasstransportBoundary
ElementMatrix * CreateKMatrixMasstransportBoundary(void)
Definition: Numericalflux.cpp:676
Hook::configure
void configure(DataSet *dataset)
Definition: Hook.cpp:145
Numericalflux::CreateKMatrixAdjointBalancethickness
ElementMatrix * CreateKMatrixAdjointBalancethickness(void)
Definition: Numericalflux.cpp:467
Numericalflux::nodes
Node ** nodes
Definition: Numericalflux.h:33
GsetEnum
@ GsetEnum
Definition: EnumDefinitions.h:1093
Numericalflux::CreateKMatrixMasstransport
ElementMatrix * CreateKMatrixMasstransport(void)
Definition: Numericalflux.cpp:664
UNDEF
#define UNDEF
Definition: constants.h:8
MARSHALLING
#define MARSHALLING(FIELD)
Definition: Marshalling.h:29
Numericalflux::GetNumberOfNodesOneSide
int GetNumberOfNodesOneSide(void)
Definition: Numericalflux.cpp:360
GaussTria::GaussPoint
void GaussPoint(int ig)
Definition: GaussTria.cpp:479
ElementMatrix::Transpose
void Transpose(void)
Definition: ElementMatrix.cpp:502
GetVerticesCoordinates
void GetVerticesCoordinates(IssmDouble *xyz, Vertex **vertices, int numvertices, bool spherical)
Definition: Vertex.cpp:225
Numericalflux::CreatePVectorBalancethickness
ElementVector * CreatePVectorBalancethickness(void)
Definition: Numericalflux.cpp:843
Tria::FiniteElement
int FiniteElement(void)
Definition: Tria.cpp:1318
Hook::copy
Object * copy(void)
Definition: Hook.cpp:61
Numericalflux::CreatePVectorBalancethicknessBoundary
ElementVector * CreatePVectorBalancethicknessBoundary(void)
Definition: Numericalflux.cpp:855
Input2
Definition: Input2.h:18
MARSHALLING_BACKWARD
@ MARSHALLING_BACKWARD
Definition: Marshalling.h:10
Numericalflux::hvertices
Hook * hvertices
Definition: Numericalflux.h:28
Numericalflux::CreatePVectorMasstransport
ElementVector * CreatePVectorMasstransport(void)
Definition: Numericalflux.cpp:923
Element::IsIceInElement
bool IsIceInElement()
Definition: Element.cpp:2021
Node
Definition: Node.h:23
Node::Lid
int Lid(void)
Definition: Node.cpp:618
Numericalflux
Definition: Numericalflux.h:18
_error_
#define _error_(StreamArgs)
Definition: exceptions.h:49
Tria
Definition: Tria.h:29
InternalEnum
@ InternalEnum
Definition: EnumDefinitions.h:1131
IoModel::faces
int * faces
Definition: IoModel.h:88
Numericalflux::CreatePVectorBalancethicknessInternal
ElementVector * CreatePVectorBalancethicknessInternal(void)
Definition: Numericalflux.cpp:916
VxEnum
@ VxEnum
Definition: EnumDefinitions.h:846
Numericalflux::vertices
Vertex ** vertices
Definition: Numericalflux.h:32
Numericalflux::parameters
Parameters * parameters
Definition: Numericalflux.h:34
GaussTria::GaussEdgeCenter
void GaussEdgeCenter(int index1, int index2)
Definition: GaussTria.cpp:423
Parameters::FindParam
void FindParam(bool *pinteger, int enum_type)
Definition: Parameters.cpp:262
ThicknessEnum
@ ThicknessEnum
Definition: EnumDefinitions.h:840
ElementVector::AddToGlobal
void AddToGlobal(Vector< IssmDouble > *pf)
Definition: ElementVector.cpp:155
AnalysisTypeEnum
@ AnalysisTypeEnum
Definition: EnumDefinitions.h:36
ElementVector
Definition: ElementVector.h:20
ElementMatrix::AddToGlobal
void AddToGlobal(Matrix< IssmDouble > *Kff, Matrix< IssmDouble > *Kfs)
Definition: ElementMatrix.cpp:271
IoModel::elements
int * elements
Definition: IoModel.h:79
Hook::Echo
void Echo(void)
Definition: Hook.cpp:104
Input2::GetInputValue
virtual void GetInputValue(IssmDouble *pvalue, Gauss *gauss)
Definition: Input2.h:38
Vertex
Definition: Vertex.h:19
Gauss::weight
IssmDouble weight
Definition: Gauss.h:11
FsetEnum
@ FsetEnum
Definition: EnumDefinitions.h:1075
DataSet
Declaration of DataSet class.
Definition: DataSet.h:14
ElementMatrix
Definition: ElementMatrix.h:19
Numericalflux::CreateKMatrixBalancethicknessBoundary
ElementMatrix * CreateKMatrixBalancethicknessBoundary(void)
Definition: Numericalflux.cpp:505
Tria::GetInput2
Input2 * GetInput2(int enumtype)
Definition: Tria.cpp:1880
Tria::GetVertexIndex
int GetVertexIndex(Vertex *vertex)
Definition: Tria.cpp:2368
Hook::Marshall
void Marshall(char **pmarshalled_data, int *pmarshalled_data_size, int marshall_direction)
Definition: Hook.cpp:122
ElementMatrix::values
IssmDouble * values
Definition: ElementMatrix.h:26
Generated on Thu Jul 2 2020 08:09:29 for Ice Sheet System Model by   doxygen 1.8.19