source: issm/trunk/src/c/classes/Loads/Numericalflux.cpp@ 24313

Last change on this file since 24313 was 24313, checked in by Mathieu Morlighem, 5 years ago

merged trunk-jpl and trunk for revision 24310
File size: 30.6 KB
Line 
1/*!\file Numericalflux.c
2 * \brief: implementation of the Numericalflux object
3 */
4
5/*Headers:*/
6/*{{{*/
7#ifdef HAVE_CONFIG_H
8#include <config.h>
9#else
10#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!"
11#endif
12
13#include "shared/shared.h"
14#include "../classes.h"
15/*}}}*/
16
17/*Load macros*/
18#define NUMVERTICES 2
19
20/*Numericalflux constructors and destructor*/
21Numericalflux::Numericalflux(){/*{{{*/
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}
29/*}}}*/
30Numericalflux::Numericalflux(int numericalflux_id,int i,int index,IoModel* iomodel){/*{{{*/
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}
135/*}}}*/
136Numericalflux::~Numericalflux(){/*{{{*/
137 this->parameters=NULL;
138 delete helement;
139 delete hnodes;
140 delete hvertices;
141}
142/*}}}*/
143
144/*Object virtual functions definitions:*/
145Object* Numericalflux::copy() {/*{{{*/
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}
171/*}}}*/
172void Numericalflux::DeepEcho(void){/*{{{*/
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}
187/*}}}*/
188void Numericalflux::Echo(void){/*{{{*/
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}
198/*}}}*/
199int Numericalflux::Id(void){/*{{{*/
200 return id;
201}
202/*}}}*/
203void Numericalflux::Marshall(char** pmarshalled_data,int* pmarshalled_data_size, int marshall_direction){ /*{{{*/
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}
229/*}}}*/
230int Numericalflux::ObjectEnum(void){/*{{{*/
231
232 return NumericalfluxEnum;
233
234}
235/*}}}*/
236
237/*Load virtual functions definitions:*/
238void Numericalflux::Configure(Elements* elementsin,Loads* loadsin,Nodes* nodesin,Vertices* verticesin,Materials* materialsin,Parameters* parametersin){/*{{{*/
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}
254/*}}}*/
255void Numericalflux::CreateKMatrix(Matrix<IssmDouble>* Kff, Matrix<IssmDouble>* Kfs){/*{{{*/
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}
284/*}}}*/
285void Numericalflux::CreatePVector(Vector<IssmDouble>* pf){/*{{{*/
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}
313/*}}}*/
314void Numericalflux::GetNodesLidList(int* lidlist){/*{{{*/
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}
322/*}}}*/
323void Numericalflux::GetNodesSidList(int* sidlist){/*{{{*/
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}
331/*}}}*/
332int Numericalflux::GetNumberOfNodes(void){/*{{{*/
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}
359/*}}}*/
360int Numericalflux::GetNumberOfNodesOneSide(void){/*{{{*/
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}
373/*}}}*/
374bool Numericalflux::IsPenalty(void){/*{{{*/
375 return false;
376}
377/*}}}*/
378void Numericalflux::PenaltyCreateKMatrix(Matrix<IssmDouble>* Kff, Matrix<IssmDouble>* Kfs,IssmDouble kmax){/*{{{*/
379
380 /*No stiffness loads applied, do nothing: */
381 return;
382
383}
384/*}}}*/
385void Numericalflux::PenaltyCreatePVector(Vector<IssmDouble>* pf,IssmDouble kmax){/*{{{*/
386
387 /*No penalty loads applied, do nothing: */
388 return;
389
390}
391/*}}}*/
392void Numericalflux::ResetHooks(){/*{{{*/
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}
405/*}}}*/
406void Numericalflux::SetCurrentConfiguration(Elements* elementsin,Loads* loadsin,Nodes* nodesin,Vertices* verticesin,Materials* materialsin,Parameters* parametersin){/*{{{*/
407 /*Nothing to do :)*/
408
409}
410/*}}}*/
411void Numericalflux::SetwiseNodeConnectivity(int* pd_nz,int* po_nz,Node* node,bool* flags,int* flagsindices,int set1_enum,int set2_enum){/*{{{*/
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}
464/*}}}*/
465
466/*Numericalflux management*/
467ElementMatrix* Numericalflux::CreateKMatrixAdjointBalancethickness(void){/*{{{*/
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}
478/*}}}*/
479ElementMatrix* Numericalflux::CreateKMatrixAdjointBalancethicknessBoundary(void){/*{{{*/
480
481 ElementMatrix* Ke=CreateKMatrixBalancethicknessBoundary();
482 if(Ke) Ke->Transpose();
483 return Ke;
484}
485/*}}}*/
486ElementMatrix* Numericalflux::CreateKMatrixAdjointBalancethicknessInternal(void){/*{{{*/
487
488 ElementMatrix* Ke=CreateKMatrixBalancethicknessInternal();
489 if (Ke) Ke->Transpose();
490 return Ke;
491}
492/*}}}*/
493ElementMatrix* Numericalflux::CreateKMatrixBalancethickness(void){/*{{{*/
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}
504/*}}}*/
505ElementMatrix* Numericalflux::CreateKMatrixBalancethicknessBoundary(void){/*{{{*/
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 Input* vxaverage_input=tria->inputs->GetInput(VxEnum);
519 Input* vyaverage_input=tria->inputs->GetInput(VyEnum);
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}
568/*}}}*/
569ElementMatrix* Numericalflux::CreateKMatrixBalancethicknessInternal(void){/*{{{*/
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 Input* vxaverage_input=tria->inputs->GetInput(VxEnum);
594 Input* vyaverage_input=tria->inputs->GetInput(VyEnum);
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}
663/*}}}*/
664ElementMatrix* Numericalflux::CreateKMatrixMasstransport(void){/*{{{*/
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}
675/*}}}*/
676ElementMatrix* Numericalflux::CreateKMatrixMasstransportBoundary(void){/*{{{*/
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 Input* vxaverage_input=tria->inputs->GetInput(VxEnum); _assert_(vxaverage_input);
691 Input* vyaverage_input=tria->inputs->GetInput(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}
740/*}}}*/
741ElementMatrix* Numericalflux::CreateKMatrixMasstransportInternal(void){/*{{{*/
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 Input* vxaverage_input=tria->inputs->GetInput(VxEnum);
767 Input* vyaverage_input=tria->inputs->GetInput(VyEnum);
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}
836/*}}}*/
837ElementVector* Numericalflux::CreatePVectorAdjointBalancethickness(void){/*{{{*/
838
839 /*No PVector for the Adjoint*/
840 return NULL;
841}
842/*}}}*/
843ElementVector* Numericalflux::CreatePVectorBalancethickness(void){/*{{{*/
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}
854/*}}}*/
855ElementVector* Numericalflux::CreatePVectorBalancethicknessBoundary(void){/*{{{*/
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 Input* vxaverage_input = tria->inputs->GetInput(VxEnum); _assert_(vxaverage_input);
869 Input* vyaverage_input = tria->inputs->GetInput(VyEnum); _assert_(vyaverage_input);
870 Input* thickness_input = tria->inputs->GetInput(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}
915/*}}}*/
916ElementVector* Numericalflux::CreatePVectorBalancethicknessInternal(void){/*{{{*/
917
918 /*Nothing added to PVector*/
919 return NULL;
920
921}
922/*}}}*/
923ElementVector* Numericalflux::CreatePVectorMasstransport(void){/*{{{*/
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}
934/*}}}*/
935ElementVector* Numericalflux::CreatePVectorMasstransportBoundary(void){/*{{{*/
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 Input* vxaverage_input = tria->inputs->GetInput(VxEnum); _assert_(vxaverage_input);
950 Input* vyaverage_input = tria->inputs->GetInput(VyEnum); _assert_(vyaverage_input);
951 Input* spcthickness_input = tria->inputs->GetInput(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}
997/*}}}*/
998ElementVector* Numericalflux::CreatePVectorMasstransportInternal(void){/*{{{*/
999
1000 /*Nothing added to PVector*/
1001 return NULL;
1002
1003}
1004/*}}}*/
1005void Numericalflux::GetNormal(IssmDouble* normal,IssmDouble xyz_list[4][3]){/*{{{*/
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}
1018/*}}}*/
Note: See TracBrowser for help on using the repository browser.