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

Last change on this file since 25836 was 25836, checked in by Mathieu Morlighem, 4 years ago

merged trunk-jpl and trunk for revision 25834
File size: 30.2 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(MarshallHandle* marshallhandle){ /*{{{*/
204
205 _assert_(this);
206
207 /*ok, marshall operations: */
208 int object_enum = NumericalfluxEnum;
209 marshallhandle->call(object_enum);
210 marshallhandle->call(this->id);
211 marshallhandle->call(this->flux_type);
212 marshallhandle->call(this->flux_degree);
213
214 if(marshallhandle->OperationNumber()==MARSHALLING_LOAD){
215 this->hnodes = new Hook();
216 this->hvertices = new Hook();
217 this->helement = new Hook();
218 }
219
220 this->hnodes->Marshall(marshallhandle);
221 this->helement->Marshall(marshallhandle);
222 this->hvertices->Marshall(marshallhandle);
223
224 /*corresponding fields*/
225 nodes =(Node**)this->hnodes->deliverp();
226 vertices =(Vertex**)this->hvertices->deliverp();
227 element =(Element*)this->helement->delivers();
228
229}
230/*}}}*/
231int Numericalflux::ObjectEnum(void){/*{{{*/
232
233 return NumericalfluxEnum;
234
235}
236/*}}}*/
237
238/*Load virtual functions definitions:*/
239void Numericalflux::Configure(Elements* elementsin,Loads* loadsin,Nodes* nodesin,Vertices* verticesin,Materials* materialsin,Parameters* parametersin){/*{{{*/
240
241 /*Take care of hooking up all objects for this element, ie links the objects in the hooks to their respective
242 * datasets, using internal ids and offsets hidden in hooks: */
243 hnodes->configure((DataSet*)nodesin);
244 hvertices->configure((DataSet*)verticesin);
245 helement->configure((DataSet*)elementsin);
246
247 /*Initialize hooked fields*/
248 this->nodes = (Node**)hnodes->deliverp();
249 this->vertices = (Vertex**)hvertices->deliverp();
250 this->element = (Element*)helement->delivers();
251
252 /*point parameters to real dataset: */
253 this->parameters=parametersin;
254}
255/*}}}*/
256void Numericalflux::CreateKMatrix(Matrix<IssmDouble>* Kff, Matrix<IssmDouble>* Kfs){/*{{{*/
257
258 /*recover some parameters*/
259 ElementMatrix* Ke=NULL;
260 int analysis_type;
261 this->parameters->FindParam(&analysis_type,AnalysisTypeEnum);
262
263 /*Just branch to the correct element stiffness matrix generator, according to the type of analysis we are carrying out: */
264 switch(analysis_type){
265 case MasstransportAnalysisEnum:
266 Ke=CreateKMatrixMasstransport();
267 break;
268 case BalancethicknessAnalysisEnum:
269 Ke=CreateKMatrixBalancethickness();
270 break;
271 case AdjointBalancethicknessAnalysisEnum:
272 Ke=CreateKMatrixAdjointBalancethickness();
273 break;
274 default:
275 _error_("analysis " << analysis_type << " (" << EnumToStringx(analysis_type) << ") not supported yet");
276 }
277
278 /*Add to global matrix*/
279 if(Ke){
280 Ke->AddToGlobal(Kff,Kfs);
281 delete Ke;
282 }
283
284}
285/*}}}*/
286void Numericalflux::CreatePVector(Vector<IssmDouble>* pf){/*{{{*/
287
288 /*recover some parameters*/
289 ElementVector* pe=NULL;
290 int analysis_type;
291 this->parameters->FindParam(&analysis_type,AnalysisTypeEnum);
292
293 switch(analysis_type){
294 case MasstransportAnalysisEnum:
295 pe=CreatePVectorMasstransport();
296 break;
297 case BalancethicknessAnalysisEnum:
298 pe=CreatePVectorBalancethickness();
299 break;
300 case AdjointBalancethicknessAnalysisEnum:
301 pe=CreatePVectorAdjointBalancethickness();
302 break;
303 default:
304 _error_("analysis " << analysis_type << " (" << EnumToStringx(analysis_type) << ") not supported yet");
305 }
306
307 /*Add to global matrix*/
308 if(pe){
309 pe->AddToGlobal(pf);
310 delete pe;
311 }
312
313}
314/*}}}*/
315void Numericalflux::GetNodesLidList(int* lidlist){/*{{{*/
316
317 _assert_(lidlist);
318 _assert_(nodes);
319
320 int numnodes = this->GetNumberOfNodes();
321 for(int i=0;i<numnodes;i++) lidlist[i]=nodes[i]->Lid();
322}
323/*}}}*/
324void Numericalflux::GetNodesSidList(int* sidlist){/*{{{*/
325
326 _assert_(sidlist);
327 _assert_(nodes);
328
329 int numnodes = this->GetNumberOfNodes();
330 for(int i=0;i<numnodes;i++) sidlist[i]=nodes[i]->Sid();
331}
332/*}}}*/
333int Numericalflux::GetNumberOfNodes(void){/*{{{*/
334
335 if(this->flux_degree==P0DGEnum){
336 switch(this->flux_type){
337 case InternalEnum:
338 return 2;
339 case BoundaryEnum:
340 return 1;
341 default:
342 _error_("Numericalflux type " << EnumToStringx(this->flux_type) << " not supported yet");
343 }
344 }
345 else if(this->flux_degree==P1DGEnum){
346 switch(this->flux_type){
347 case InternalEnum:
348 return 4;
349 case BoundaryEnum:
350 return 2;
351 default:
352 _error_("Numericalflux type " << EnumToStringx(this->flux_type) << " not supported yet");
353 }
354 }
355 else{
356 _error_("Numericalflux " << EnumToStringx(this->flux_degree) << " not supported yet");
357 }
358
359}
360/*}}}*/
361int Numericalflux::GetNumberOfNodesOneSide(void){/*{{{*/
362
363 if(this->flux_degree==P0DGEnum){
364 return 1;
365 }
366 else if(this->flux_degree==P1DGEnum){
367 return 2;
368 }
369 else{
370 _error_("Numericalflux " << EnumToStringx(this->flux_degree) << " not supported yet");
371 }
372
373}
374/*}}}*/
375bool Numericalflux::IsPenalty(void){/*{{{*/
376 return false;
377}
378/*}}}*/
379void Numericalflux::PenaltyCreateKMatrix(Matrix<IssmDouble>* Kff, Matrix<IssmDouble>* Kfs,IssmDouble kmax){/*{{{*/
380
381 /*No stiffness loads applied, do nothing: */
382 return;
383
384}
385/*}}}*/
386void Numericalflux::PenaltyCreatePVector(Vector<IssmDouble>* pf,IssmDouble kmax){/*{{{*/
387
388 /*No penalty loads applied, do nothing: */
389 return;
390
391}
392/*}}}*/
393void Numericalflux::ResetHooks(){/*{{{*/
394
395 this->nodes = NULL;
396 this->vertices = NULL;
397 this->element = NULL;
398 this->parameters = NULL;
399
400 /*Get Element type*/
401 this->hnodes->reset();
402 this->hvertices->reset();
403 this->helement->reset();
404
405}
406/*}}}*/
407void Numericalflux::SetCurrentConfiguration(Elements* elementsin,Loads* loadsin,Nodes* nodesin,Vertices* verticesin,Materials* materialsin,Parameters* parametersin){/*{{{*/
408 /*Nothing to do :)*/
409
410}
411/*}}}*/
412void Numericalflux::SetwiseNodeConnectivity(int* pd_nz,int* po_nz,Node* node,bool* flags,int* flagsindices,int* flagsindices_counter,int set1_enum,int set2_enum){/*{{{*/
413
414 /*Output */
415 int d_nz = 0;
416 int o_nz = 0;
417
418 /*Loop over all nodes*/
419 for(int i=0;i<this->GetNumberOfNodes();i++){
420
421 if(!flags[this->nodes[i]->Lid()]){
422
423 /*flag current node so that no other element processes it*/
424 flags[this->nodes[i]->Lid()]=true;
425
426 flagsindices[flagsindices_counter[0]]=this->nodes[i]->Lid();
427 flagsindices_counter[0]++;
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->GetInput(VxEnum); _assert_(vxaverage_input);
519 Input* vyaverage_input=tria->GetInput(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 while(gauss->next()){
545
546 tria->GetSegmentNodalFunctions(&basis[0],gauss,index1,index2,tria->FiniteElement());
547
548 vxaverage_input->GetInputValue(&vx,gauss);
549 vyaverage_input->GetInputValue(&vy,gauss);
550 UdotN=vx*normal[0]+vy*normal[1];
551 tria->GetSegmentJacobianDeterminant(&Jdet,&xyz_list[0][0],gauss);
552 DL=gauss->weight*Jdet*UdotN;
553
554 for(int i=0;i<numnodes;i++){
555 for(int j=0;j<numnodes;j++){
556 Ke->values[i*numnodes+j]+=DL*basis[i]*basis[j];
557 }
558 }
559 }
560
561 /*Clean up and return*/
562 xDelete<IssmDouble>(basis);
563 delete gauss;
564 return Ke;
565}
566/*}}}*/
567ElementMatrix* Numericalflux::CreateKMatrixBalancethicknessInternal(void){/*{{{*/
568
569 /*Initialize Element matrix and return if necessary*/
570 Tria* tria=(Tria*)element;
571 if(!tria->IsIceInElement()) return NULL;
572
573 /* Intermediaries*/
574 IssmDouble A1,A2,Jdet,vx,vy,UdotN;
575 IssmDouble xyz_list[NUMVERTICES][3];
576 IssmDouble normal[2];
577
578 /*Fetch number of nodes for this flux*/
579 int numnodes = this->GetNumberOfNodes();
580 int numnodes_plus = this->GetNumberOfNodesOneSide();
581 int numnodes_minus = numnodes_plus; /*For now we are not doing p-adaptive DG*/
582 _assert_(numnodes==numnodes_plus+numnodes_minus);
583
584 /*Initialize variables*/
585 ElementMatrix *Ke = new ElementMatrix(nodes,numnodes,this->parameters);
586 IssmDouble *basis_plus = xNew<IssmDouble>(numnodes_plus);
587 IssmDouble *basis_minus = xNew<IssmDouble>(numnodes_minus);
588
589 /*Retrieve all inputs and parameters*/
590 GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
591 Input* vxaverage_input=tria->GetInput(VxEnum); _assert_(vxaverage_input);
592 Input* vyaverage_input=tria->GetInput(VyEnum); _assert_(vyaverage_input);
593 GetNormal(&normal[0],xyz_list);
594
595 /* Start looping on the number of gaussian points: */
596 int index1=tria->GetVertexIndex(vertices[0]);
597 int index2=tria->GetVertexIndex(vertices[1]);
598 GaussTria* gauss=new GaussTria(index1,index2,2);
599 while(gauss->next()){
600
601 tria->GetSegmentNodalFunctions(&basis_plus[0] ,gauss,index1,index2,tria->FiniteElement());
602 tria->GetSegmentNodalFunctions(&basis_minus[0],gauss,index1,index2,tria->FiniteElement());
603
604 vxaverage_input->GetInputValue(&vx,gauss);
605 vyaverage_input->GetInputValue(&vy,gauss);
606 UdotN=vx*normal[0]+vy*normal[1];
607 tria->GetSegmentJacobianDeterminant(&Jdet,&xyz_list[0][0],gauss);
608 A1=gauss->weight*Jdet*UdotN/2;
609 A2=gauss->weight*Jdet*fabs(UdotN)/2;
610
611 /*Term 1 (numerical flux): {Hv}.[[phi]] = 0.5(H+v+ + H-v-)(phi+n+ + phi-n-)
612 * = v.n/2 (H+phi+ + H-phi+ -H+phi- -H-phi-)
613 * = v.n/2 (H+phi+ + H-phi+ -H+phi- -H-phi-)
614 *
615 *Term 2 (stabilization) |v.n|/2 [[H]].[[phi]] = |v.n|/2 (H+n+ + H-n-)(phi+n+ + phi-n-)
616 * = |v.n|/2 (H+phi+ -H-phi+ -H+phi- +H-phi-)
617 * | A++ | A+- |
618 * K = |-----------|
619 * | A-+ | A-- |
620 *
621 *These 4 terms for each expressions are added independently*/
622
623 /*First term A++*/
624 for(int i=0;i<numnodes_plus;i++){
625 for(int j=0;j<numnodes_plus;j++){
626 Ke->values[i*numnodes+j] += A1*(basis_plus[j]*basis_plus[i]);
627 Ke->values[i*numnodes+j] += A2*(basis_plus[j]*basis_plus[i]);
628 }
629 }
630 /*Second term A+-*/
631 for(int i=0;i<numnodes_plus;i++){
632 for(int j=0;j<numnodes_minus;j++){
633 Ke->values[i*numnodes+numnodes_plus+j] += A1*(basis_minus[j]*basis_plus[i]);
634 Ke->values[i*numnodes+numnodes_plus+j] += -A2*(basis_minus[j]*basis_plus[i]);
635 }
636 }
637 /*Third term A-+*/
638 for(int i=0;i<numnodes_minus;i++){
639 for(int j=0;j<numnodes_plus;j++){
640 Ke->values[(numnodes_plus+i)*numnodes+j] += -A1*(basis_plus[j]*basis_minus[i]);
641 Ke->values[(numnodes_plus+i)*numnodes+j] += -A2*(basis_plus[j]*basis_minus[i]);
642 }
643 }
644 /*Fourth term A-+*/
645 for(int i=0;i<numnodes_minus;i++){
646 for(int j=0;j<numnodes_minus;j++){
647 Ke->values[(numnodes_plus+i)*numnodes+numnodes_plus+j] += -A1*(basis_minus[j]*basis_minus[i]);
648 Ke->values[(numnodes_plus+i)*numnodes+numnodes_plus+j] += A2*(basis_minus[j]*basis_minus[i]);
649 }
650 }
651 }
652
653 /*Clean up and return*/
654 xDelete<IssmDouble>(basis_plus);
655 xDelete<IssmDouble>(basis_minus);
656 delete gauss;
657 return Ke;
658}
659/*}}}*/
660ElementMatrix* Numericalflux::CreateKMatrixMasstransport(void){/*{{{*/
661
662 switch(this->flux_type){
663 case InternalEnum:
664 return CreateKMatrixMasstransportInternal();
665 case BoundaryEnum:
666 return CreateKMatrixMasstransportBoundary();
667 default:
668 _error_("type not supported yet");
669 }
670}
671/*}}}*/
672ElementMatrix* Numericalflux::CreateKMatrixMasstransportBoundary(void){/*{{{*/
673
674 /*Initialize Element matrix and return if necessary*/
675 Tria* tria=(Tria*)element;
676 if(!tria->IsIceInElement()) return NULL;
677
678 /* Intermediaries*/
679 IssmDouble DL,Jdet,vx,vy,mean_vx,mean_vy;
680 IssmDouble xyz_list[NUMVERTICES][3];
681 IssmDouble normal[2];
682
683 /*Retrieve all inputs and parameters*/
684 GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
685 IssmDouble dt = parameters->FindParam(TimesteppingTimeStepEnum);
686 Input* vxaverage_input=tria->GetInput(VxEnum); _assert_(vxaverage_input);
687 Input* vyaverage_input=tria->GetInput(VyEnum); _assert_(vyaverage_input);
688 GetNormal(&normal[0],xyz_list);
689
690 /*Check wether it is an inflow or outflow BC (0 is the middle of the segment)*/
691 int index1=tria->GetVertexIndex(vertices[0]);
692 int index2=tria->GetVertexIndex(vertices[1]);
693
694 GaussTria* gauss=new GaussTria();
695 gauss->GaussEdgeCenter(index1,index2);
696 vxaverage_input->GetInputValue(&mean_vx,gauss);
697 vyaverage_input->GetInputValue(&mean_vy,gauss);
698 delete gauss;
699
700 IssmDouble UdotN=mean_vx*normal[0]+mean_vy*normal[1];
701 if(UdotN<=0){
702 return NULL; /*(u,n)<0 -> inflow, PenaltyCreatePVector will take care of it*/
703 }
704
705 /*Initialize Element vector and other vectors*/
706 int numnodes = this->GetNumberOfNodes();
707 ElementMatrix *Ke = new ElementMatrix(nodes,numnodes,this->parameters);
708 IssmDouble *basis = xNew<IssmDouble>(numnodes);
709
710 /* Start looping on the number of gaussian points: */
711 gauss=new GaussTria(index1,index2,2);
712 while(gauss->next()){
713
714 tria->GetSegmentNodalFunctions(&basis[0],gauss,index1,index2,tria->FiniteElement());
715
716 vxaverage_input->GetInputValue(&vx,gauss);
717 vyaverage_input->GetInputValue(&vy,gauss);
718 UdotN=vx*normal[0]+vy*normal[1];
719 tria->GetSegmentJacobianDeterminant(&Jdet,&xyz_list[0][0],gauss);
720 DL=gauss->weight*Jdet*dt*UdotN;
721
722 for(int i=0;i<numnodes;i++){
723 for(int j=0;j<numnodes;j++){
724 Ke->values[i*numnodes+j]+=DL*basis[i]*basis[j];
725 }
726 }
727 }
728
729 /*Clean up and return*/
730 xDelete<IssmDouble>(basis);
731 delete gauss;
732 return Ke;
733}
734/*}}}*/
735ElementMatrix* Numericalflux::CreateKMatrixMasstransportInternal(void){/*{{{*/
736
737 /*Initialize Element matrix and return if necessary*/
738 Tria* tria=(Tria*)element;
739 if(!tria->IsIceInElement()) return NULL;
740
741 /* Intermediaries*/
742 IssmDouble A1,A2,Jdet,vx,vy,UdotN;
743 IssmDouble xyz_list[NUMVERTICES][3];
744 IssmDouble normal[2];
745
746 /*Fetch number of nodes for this flux*/
747 int numnodes = this->GetNumberOfNodes();
748 int numnodes_plus = this->GetNumberOfNodesOneSide();
749 int numnodes_minus = numnodes_plus; /*For now we are not doing p-adaptive DG*/
750 _assert_(numnodes==numnodes_plus+numnodes_minus);
751
752 /*Initialize variables*/
753 ElementMatrix *Ke = new ElementMatrix(nodes,numnodes,this->parameters);
754 IssmDouble *basis_plus = xNew<IssmDouble>(numnodes_plus);
755 IssmDouble *basis_minus = xNew<IssmDouble>(numnodes_minus);
756
757 /*Retrieve all inputs and parameters*/
758 GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
759 IssmDouble dt = parameters->FindParam(TimesteppingTimeStepEnum);
760 Input* vxaverage_input=tria->GetInput(VxEnum); _assert_(vxaverage_input);
761 Input* vyaverage_input=tria->GetInput(VyEnum); _assert_(vyaverage_input);
762 GetNormal(&normal[0],xyz_list);
763
764 /* Start looping on the number of gaussian points: */
765 int index1=tria->GetVertexIndex(vertices[0]);
766 int index2=tria->GetVertexIndex(vertices[1]);
767 GaussTria* gauss=new GaussTria(index1,index2,2);
768 while(gauss->next()){
769
770 tria->GetSegmentNodalFunctions(&basis_plus[0] ,gauss,index1,index2,tria->FiniteElement());
771 tria->GetSegmentNodalFunctions(&basis_minus[0],gauss,index1,index2,tria->FiniteElement());
772
773 vxaverage_input->GetInputValue(&vx,gauss);
774 vyaverage_input->GetInputValue(&vy,gauss);
775 UdotN=vx*normal[0]+vy*normal[1];
776 tria->GetSegmentJacobianDeterminant(&Jdet,&xyz_list[0][0],gauss);
777 A1=gauss->weight*Jdet*dt*UdotN/2;
778 A2=gauss->weight*Jdet*dt*fabs(UdotN)/2;
779
780 /*Term 1 (numerical flux): {Hv}.[[phi]] = 0.5(H+v+ + H-v-)(phi+n+ + phi-n-)
781 * = v.n/2 (H+phi+ + H-phi+ -H+phi- -H-phi-)
782 * = v.n/2 (H+phi+ + H-phi+ -H+phi- -H-phi-)
783 *
784 *Term 2 (stabilization) |v.n|/2 [[H]].[[phi]] = |v.n|/2 (H+n+ + H-n-)(phi+n+ + phi-n-)
785 * = |v.n|/2 (H+phi+ -H-phi+ -H+phi- +H-phi-)
786 * | A++ | A+- |
787 * K = |-----------|
788 * | A-+ | A-- |
789 *
790 *These 4 terms for each expressions are added independently*/
791
792 /*First term A++*/
793 for(int i=0;i<numnodes_plus;i++){
794 for(int j=0;j<numnodes_plus;j++){
795 Ke->values[i*numnodes+j] += A1*(basis_plus[j]*basis_plus[i]);
796 Ke->values[i*numnodes+j] += A2*(basis_plus[j]*basis_plus[i]);
797 }
798 }
799 /*Second term A+-*/
800 for(int i=0;i<numnodes_plus;i++){
801 for(int j=0;j<numnodes_minus;j++){
802 Ke->values[i*numnodes+numnodes_plus+j] += A1*(basis_minus[j]*basis_plus[i]);
803 Ke->values[i*numnodes+numnodes_plus+j] += -A2*(basis_minus[j]*basis_plus[i]);
804 }
805 }
806 /*Third term A-+*/
807 for(int i=0;i<numnodes_minus;i++){
808 for(int j=0;j<numnodes_plus;j++){
809 Ke->values[(numnodes_plus+i)*numnodes+j] += -A1*(basis_plus[j]*basis_minus[i]);
810 Ke->values[(numnodes_plus+i)*numnodes+j] += -A2*(basis_plus[j]*basis_minus[i]);
811 }
812 }
813 /*Fourth term A-+*/
814 for(int i=0;i<numnodes_minus;i++){
815 for(int j=0;j<numnodes_minus;j++){
816 Ke->values[(numnodes_plus+i)*numnodes+numnodes_plus+j] += -A1*(basis_minus[j]*basis_minus[i]);
817 Ke->values[(numnodes_plus+i)*numnodes+numnodes_plus+j] += A2*(basis_minus[j]*basis_minus[i]);
818 }
819 }
820 }
821
822 /*Clean up and return*/
823 xDelete<IssmDouble>(basis_plus);
824 xDelete<IssmDouble>(basis_minus);
825 delete gauss;
826 return Ke;
827}
828/*}}}*/
829ElementVector* Numericalflux::CreatePVectorAdjointBalancethickness(void){/*{{{*/
830
831 /*No PVector for the Adjoint*/
832 return NULL;
833}
834/*}}}*/
835ElementVector* Numericalflux::CreatePVectorBalancethickness(void){/*{{{*/
836
837 switch(this->flux_type){
838 case InternalEnum:
839 return CreatePVectorBalancethicknessInternal();
840 case BoundaryEnum:
841 return CreatePVectorBalancethicknessBoundary();
842 default:
843 _error_("type not supported yet");
844 }
845}
846/*}}}*/
847ElementVector* Numericalflux::CreatePVectorBalancethicknessBoundary(void){/*{{{*/
848
849 /*Initialize Load Vector and return if necessary*/
850 Tria* tria=(Tria*)element;
851 if(!tria->IsIceInElement()) return NULL;
852
853 /* Intermediaries*/
854 IssmDouble DL,Jdet,vx,vy,mean_vx,mean_vy,thickness;
855 IssmDouble xyz_list[NUMVERTICES][3];
856 IssmDouble normal[2];
857
858 /*Retrieve all inputs and parameters*/
859 GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
860 Input* vxaverage_input = tria->GetInput(VxEnum); _assert_(vxaverage_input);
861 Input* vyaverage_input = tria->GetInput(VyEnum); _assert_(vyaverage_input);
862 Input* thickness_input = tria->GetInput(ThicknessEnum); _assert_(thickness_input);
863 GetNormal(&normal[0],xyz_list);
864
865 /*Check wether it is an inflow or outflow BC (0 is the middle of the segment)*/
866 int index1=tria->GetVertexIndex(vertices[0]);
867 int index2=tria->GetVertexIndex(vertices[1]);
868 GaussTria* gauss=new GaussTria();
869 gauss->GaussEdgeCenter(index1,index2);
870 vxaverage_input->GetInputValue(&mean_vx,gauss);
871 vyaverage_input->GetInputValue(&mean_vy,gauss);
872 delete gauss;
873 IssmDouble UdotN=mean_vx*normal[0]+mean_vy*normal[1];
874 if(UdotN>0){
875 return NULL; /*(u,n)>0 -> outflow, PenaltyCreateKMatrix will take care of it*/
876 }
877
878 /*Initialize Load Vector */
879 int numnodes = this->GetNumberOfNodes();
880 ElementVector *pe = new ElementVector(nodes,numnodes,this->parameters);
881 IssmDouble *basis = xNew<IssmDouble>(numnodes);
882
883 /* Start looping on the number of gaussian points: */
884 gauss=new GaussTria(index1,index2,2);
885 while(gauss->next()){
886
887 tria->GetSegmentNodalFunctions(&basis[0],gauss,index1,index2,tria->FiniteElement());
888
889 vxaverage_input->GetInputValue(&vx,gauss);
890 vyaverage_input->GetInputValue(&vy,gauss);
891 thickness_input->GetInputValue(&thickness,gauss);
892
893 UdotN=vx*normal[0]+vy*normal[1];
894 tria->GetSegmentJacobianDeterminant(&Jdet,&xyz_list[0][0],gauss);
895 DL= - gauss->weight*Jdet*UdotN*thickness;
896
897 for(int i=0;i<numnodes;i++) pe->values[i] += DL*basis[i];
898 }
899
900 /*Clean up and return*/
901 xDelete<IssmDouble>(basis);
902 delete gauss;
903 return pe;
904}
905/*}}}*/
906ElementVector* Numericalflux::CreatePVectorBalancethicknessInternal(void){/*{{{*/
907
908 /*Nothing added to PVector*/
909 return NULL;
910
911}
912/*}}}*/
913ElementVector* Numericalflux::CreatePVectorMasstransport(void){/*{{{*/
914
915 switch(this->flux_type){
916 case InternalEnum:
917 return CreatePVectorMasstransportInternal();
918 case BoundaryEnum:
919 return CreatePVectorMasstransportBoundary();
920 default:
921 _error_("type not supported yet");
922 }
923}
924/*}}}*/
925ElementVector* Numericalflux::CreatePVectorMasstransportBoundary(void){/*{{{*/
926
927 /*Initialize Load Vector and return if necessary*/
928 Tria* tria=(Tria*)element;
929 if(!tria->IsIceInElement()) return NULL;
930
931 /* Intermediaries*/
932 IssmDouble DL,Jdet,vx,vy,mean_vx,mean_vy,thickness;
933 IssmDouble xyz_list[NUMVERTICES][3];
934 IssmDouble normal[2];
935
936 /*Retrieve all inputs and parameters*/
937 GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
938 IssmDouble dt = parameters->FindParam(TimesteppingTimeStepEnum);
939 Input* vxaverage_input = tria->GetInput(VxEnum); _assert_(vxaverage_input);
940 Input* vyaverage_input = tria->GetInput(VyEnum); _assert_(vyaverage_input);
941 Input* spcthickness_input = tria->GetInput(MasstransportSpcthicknessEnum); _assert_(spcthickness_input);
942 GetNormal(&normal[0],xyz_list);
943
944 /*Check wether it is an inflow or outflow BC (0 is the middle of the segment)*/
945 int index1=tria->GetVertexIndex(vertices[0]);
946 int index2=tria->GetVertexIndex(vertices[1]);
947 GaussTria* gauss=new GaussTria();
948 gauss->GaussEdgeCenter(index1,index2);
949 vxaverage_input->GetInputValue(&mean_vx,gauss);
950 vyaverage_input->GetInputValue(&mean_vy,gauss);
951 delete gauss;
952 IssmDouble UdotN=mean_vx*normal[0]+mean_vy*normal[1];
953 if(UdotN>0){
954 return NULL; /*(u,n)>0 -> outflow, PenaltyCreateKMatrix will take care of it*/
955 }
956
957 /*Initialize Load Vector */
958 int numnodes = this->GetNumberOfNodes();
959 ElementVector *pe = new ElementVector(nodes,numnodes,this->parameters);
960 IssmDouble *basis = xNew<IssmDouble>(numnodes);
961
962 /* Start looping on the number of gaussian points: */
963 gauss=new GaussTria(index1,index2,2);
964 while(gauss->next()){
965
966 tria->GetSegmentNodalFunctions(&basis[0],gauss,index1,index2,tria->FiniteElement());
967
968 vxaverage_input->GetInputValue(&vx,gauss);
969 vyaverage_input->GetInputValue(&vy,gauss);
970 spcthickness_input->GetInputValue(&thickness,gauss);
971 if(xIsNan<IssmDouble>(thickness)) _error_("Cannot weakly apply constraint because NaN was provided");
972
973 UdotN=vx*normal[0]+vy*normal[1];
974 tria->GetSegmentJacobianDeterminant(&Jdet,&xyz_list[0][0],gauss);
975 DL= - gauss->weight*Jdet*dt*UdotN*thickness;
976
977 for(int i=0;i<numnodes;i++) pe->values[i] += DL*basis[i];
978 }
979
980 /*Clean up and return*/
981 xDelete<IssmDouble>(basis);
982 delete gauss;
983 return pe;
984}
985/*}}}*/
986ElementVector* Numericalflux::CreatePVectorMasstransportInternal(void){/*{{{*/
987
988 /*Nothing added to PVector*/
989 return NULL;
990
991}
992/*}}}*/
993void Numericalflux::GetNormal(IssmDouble* normal,IssmDouble xyz_list[4][3]){/*{{{*/
994
995 /*Build unit outward pointing vector*/
996 IssmDouble vector[2];
997
998 vector[0]=xyz_list[1][0] - xyz_list[0][0];
999 vector[1]=xyz_list[1][1] - xyz_list[0][1];
1000
1001 IssmDouble norm=sqrt(pow(vector[0],2.0)+pow(vector[1],2.0));
1002
1003 normal[0]= + vector[1]/norm;
1004 normal[1]= - vector[0]/norm;
1005}
1006/*}}}*/
Note: See TracBrowser for help on using the repository browser.