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

Last change on this file since 16560 was 16560, checked in by Mathieu Morlighem, 11 years ago

merged trunk-jpl and trunk for revision 16554
File size: 28.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#define NUMNODES_INTERNAL 4
20#define NUMNODES_BOUNDARY 2
21
22/*Numericalflux constructors and destructor*/
23/*FUNCTION Numericalflux::Numericalflux(){{{*/
24Numericalflux::Numericalflux(){
25 this->inputs = NULL;
26 this->parameters = NULL;
27 this->helement = NULL;
28 this->element = NULL;
29 this->hnodes = NULL;
30 this->hvertices = NULL;
31 this->nodes = NULL;
32}
33/*}}}*/
34/*FUNCTION Numericalflux::Numericalflux(int id, int i, IoModel* iomodel, int analysis_type) {{{*/
35Numericalflux::Numericalflux(int numericalflux_id,int i,int index,IoModel* iomodel, int in_analysis_type){
36
37 /* Intermediary */
38 int j;
39 int pos1,pos2,pos3,pos4;
40 int num_nodes;
41
42 /*numericalflux constructor data: */
43 int numericalflux_elem_ids[2];
44 int numericalflux_vertex_ids[2];
45 int numericalflux_node_ids[4];
46 int numericalflux_type;
47
48 /*Get edge*/
49 int i1 = iomodel->faces[4*index+0];
50 int i2 = iomodel->faces[4*index+1];
51 int e1 = iomodel->faces[4*index+2];
52 int e2 = iomodel->faces[4*index+3];
53
54 /*First, see wether this is an internal or boundary edge (if e2=-1)*/
55 if(e2==-1){
56 /* Boundary edge, only one element */
57 num_nodes=2;
58 numericalflux_type=BoundaryEnum;
59 numericalflux_elem_ids[0]=e1;
60 }
61 else{
62 /* internal edge: connected to 2 elements */
63 num_nodes=4;
64 numericalflux_type=InternalEnum;
65 numericalflux_elem_ids[0]=e1;
66 numericalflux_elem_ids[1]=e2;
67 }
68
69 /*1: Get vertices ids*/
70 numericalflux_vertex_ids[0]=i1;
71 numericalflux_vertex_ids[1]=i2;
72
73 /*2: Get node ids*/
74 if (numericalflux_type==InternalEnum){
75
76 /*Now, we must get the nodes of the 4 nodes located on the edge*/
77
78 /*2: Get the column where these ids are located in the index*/
79 pos1=pos2=pos3=pos4=UNDEF;
80 for(j=0;j<3;j++){
81 if(iomodel->elements[3*(e1-1)+j]==i1) pos1=j+1;
82 if(iomodel->elements[3*(e1-1)+j]==i2) pos2=j+1;
83 if(iomodel->elements[3*(e2-1)+j]==i1) pos3=j+1;
84 if(iomodel->elements[3*(e2-1)+j]==i2) pos4=j+1;
85 }
86 _assert_(pos1!=UNDEF && pos2!=UNDEF && pos3!=UNDEF && pos4!=UNDEF);
87
88 /*3: We have the id of the elements and the position of the vertices in the index
89 * we can compute their dofs!*/
90 numericalflux_node_ids[0]=iomodel->nodecounter+3*(e1-1)+pos1;
91 numericalflux_node_ids[1]=iomodel->nodecounter+3*(e1-1)+pos2;
92 numericalflux_node_ids[2]=iomodel->nodecounter+3*(e2-1)+pos3;
93 numericalflux_node_ids[3]=iomodel->nodecounter+3*(e2-1)+pos4;
94 }
95 else{
96
97 /*2: Get the column where these ids are located in the index*/
98 pos1=pos2=UNDEF;
99 for(j=0;j<3;j++){
100 if(iomodel->elements[3*(e1-1)+j]==i1) pos1=j+1;
101 if(iomodel->elements[3*(e1-1)+j]==i2) pos2=j+1;
102 }
103 _assert_(pos1!=UNDEF && pos2!=UNDEF);
104
105 /*3: We have the id of the elements and the position of the vertices in the index
106 * we can compute their dofs!*/
107 numericalflux_node_ids[0]=iomodel->nodecounter+3*(e1-1)+pos1;
108 numericalflux_node_ids[1]=iomodel->nodecounter+3*(e1-1)+pos2;
109 }
110
111 /*Ok, we have everything to build the object: */
112 this->id=numericalflux_id;
113 this->analysis_type=in_analysis_type;
114
115 /*Hooks: */
116 this->hnodes =new Hook(numericalflux_node_ids,num_nodes);
117 this->hvertices =new Hook(&numericalflux_vertex_ids[0],2);
118 this->helement =new Hook(numericalflux_elem_ids,1); // take only the first element for now
119
120 //intialize and add as many inputs per element as requested:
121 this->inputs=new Inputs();
122 this->inputs->AddInput(new IntInput(NumericalfluxTypeEnum,numericalflux_type));
123
124 //this->parameters: we still can't point to it, it may not even exist. Configure will handle this.
125 this->parameters=NULL;
126 this->element=NULL;
127 this->nodes=NULL;
128}
129/*}}}*/
130/*FUNCTION Numericalflux::~Numericalflux(){{{*/
131Numericalflux::~Numericalflux(){
132 delete inputs;
133 this->parameters=NULL;
134 delete helement;
135 delete hnodes;
136 delete hvertices;
137}
138/*}}}*/
139
140/*Object virtual functions definitions:*/
141/*FUNCTION Numericalflux::Echo {{{*/
142void Numericalflux::Echo(void){
143 _printf_("Numericalflux:\n");
144 _printf_(" id: " << id << "\n");
145 _printf_(" analysis_type: " << EnumToStringx(analysis_type) << "\n");
146 hnodes->Echo();
147 hvertices->Echo();
148 helement->Echo();
149 _printf_(" parameters: " << parameters << "\n");
150 _printf_(" inputs: " << inputs << "\n");
151}
152/*}}}*/
153/*FUNCTION Numericalflux::DeepEcho {{{*/
154void Numericalflux::DeepEcho(void){
155
156 _printf_("Numericalflux:\n");
157 _printf_(" id: " << id << "\n");
158 _printf_(" analysis_type: " << EnumToStringx(analysis_type) << "\n");
159 hnodes->DeepEcho();
160 hvertices->DeepEcho();
161 helement->DeepEcho();
162 _printf_(" parameters\n");
163 if(parameters)
164 parameters->DeepEcho();
165 else
166 _printf_(" NULL\n");
167 _printf_(" inputs\n");
168 inputs->DeepEcho();
169
170}
171/*}}}*/
172/*FUNCTION Numericalflux::Id {{{*/
173int Numericalflux::Id(void){
174 return id;
175}
176/*}}}*/
177/*FUNCTION Numericalflux::ObjectEnum{{{*/
178int Numericalflux::ObjectEnum(void){
179
180 return NumericalfluxEnum;
181
182}
183/*}}}*/
184/*FUNCTION Numericalflux::copy {{{*/
185Object* Numericalflux::copy() {
186
187 Numericalflux* numericalflux=NULL;
188
189 numericalflux=new Numericalflux();
190
191 /*copy fields: */
192 numericalflux->id=this->id;
193 numericalflux->analysis_type=this->analysis_type;
194 if(this->inputs){
195 numericalflux->inputs=(Inputs*)this->inputs->Copy();
196 }
197 else{
198 numericalflux->inputs=new Inputs();
199 }
200 /*point parameters: */
201 numericalflux->parameters=this->parameters;
202
203 /*now deal with hooks and objects: */
204 numericalflux->hnodes = (Hook*)this->hnodes->copy();
205 numericalflux->hvertices = (Hook*)this->hvertices->copy();
206 numericalflux->helement = (Hook*)this->helement->copy();
207
208 /*corresponding fields*/
209 numericalflux->nodes = (Node**)numericalflux->hnodes->deliverp();
210 numericalflux->vertices = (Vertex**)numericalflux->hvertices->deliverp();
211 numericalflux->element = (Element*)numericalflux->helement->delivers();
212
213 return numericalflux;
214}
215/*}}}*/
216
217/*Load virtual functions definitions:*/
218/*FUNCTION Numericalflux::Configure {{{*/
219void Numericalflux::Configure(Elements* elementsin,Loads* loadsin,Nodes* nodesin,Vertices* verticesin,Materials* materialsin,Parameters* parametersin){
220
221 /*Take care of hooking up all objects for this element, ie links the objects in the hooks to their respective
222 * datasets, using internal ids and offsets hidden in hooks: */
223 hnodes->configure((DataSet*)nodesin);
224 hvertices->configure((DataSet*)verticesin);
225 helement->configure((DataSet*)elementsin);
226
227 /*Initialize hooked fields*/
228 this->nodes = (Node**)hnodes->deliverp();
229 this->vertices = (Vertex**)hvertices->deliverp();
230 this->element = (Element*)helement->delivers();
231
232 /*point parameters to real dataset: */
233 this->parameters=parametersin;
234}
235/*}}}*/
236/*FUNCTION Numericalflux::SetCurrentConfiguration {{{*/
237void Numericalflux::SetCurrentConfiguration(Elements* elementsin,Loads* loadsin,Nodes* nodesin,Vertices* verticesin,Materials* materialsin,Parameters* parametersin){
238
239}
240/*}}}*/
241/*FUNCTION Numericalflux::CreateKMatrix {{{*/
242void Numericalflux::CreateKMatrix(Matrix<IssmDouble>* Kff, Matrix<IssmDouble>* Kfs){
243
244 /*recover some parameters*/
245 ElementMatrix* Ke=NULL;
246 int analysis_type;
247 this->parameters->FindParam(&analysis_type,AnalysisTypeEnum);
248
249 /*Just branch to the correct element stiffness matrix generator, according to the type of analysis we are carrying out: */
250 switch(analysis_type){
251 case MasstransportAnalysisEnum:
252 Ke=CreateKMatrixMasstransport();
253 break;
254 case BalancethicknessAnalysisEnum:
255 Ke=CreateKMatrixBalancethickness();
256 break;
257 case AdjointBalancethicknessAnalysisEnum:
258 Ke=CreateKMatrixAdjointBalancethickness();
259 break;
260 default:
261 _error_("analysis " << analysis_type << " (" << EnumToStringx(analysis_type) << ") not supported yet");
262 }
263
264 /*Add to global matrix*/
265 if(Ke){
266 Ke->AddToGlobal(Kff,Kfs);
267 delete Ke;
268 }
269
270}
271/*}}}*/
272/*FUNCTION Numericalflux::CreatePVector {{{*/
273void Numericalflux::CreatePVector(Vector<IssmDouble>* pf){
274
275 /*recover some parameters*/
276 ElementVector* pe=NULL;
277 int analysis_type;
278 this->parameters->FindParam(&analysis_type,AnalysisTypeEnum);
279
280 switch(analysis_type){
281 case MasstransportAnalysisEnum:
282 pe=CreatePVectorMasstransport();
283 break;
284 case BalancethicknessAnalysisEnum:
285 pe=CreatePVectorBalancethickness();
286 break;
287 case AdjointBalancethicknessAnalysisEnum:
288 pe=CreatePVectorAdjointBalancethickness();
289 break;
290 default:
291 _error_("analysis " << analysis_type << " (" << EnumToStringx(analysis_type) << ") not supported yet");
292 }
293
294 /*Add to global matrix*/
295 if(pe){
296 pe->AddToGlobal(pf);
297 delete pe;
298 }
299
300}
301/*}}}*/
302/*FUNCTION Numericalflux::GetNodesSidList{{{*/
303void Numericalflux::GetNodesSidList(int* sidlist){
304
305 int type;
306 inputs->GetInputValue(&type,NumericalfluxTypeEnum);
307 _assert_(sidlist);
308 _assert_(nodes);
309
310 switch(type){
311 case InternalEnum:
312 for(int i=0;i<NUMNODES_INTERNAL;i++) sidlist[i]=nodes[i]->Sid();
313 return;
314 case BoundaryEnum:
315 for(int i=0;i<NUMNODES_BOUNDARY;i++) sidlist[i]=nodes[i]->Sid();
316 return;
317 default:
318 _error_("Numericalflux type " << EnumToStringx(type) << " not supported yet");
319 }
320}
321/*}}}*/
322/*FUNCTION Numericalflux::GetNodesLidList{{{*/
323void Numericalflux::GetNodesLidList(int* lidlist){
324
325 int type;
326 inputs->GetInputValue(&type,NumericalfluxTypeEnum);
327 _assert_(lidlist);
328 _assert_(nodes);
329
330 switch(type){
331 case InternalEnum:
332 for(int i=0;i<NUMNODES_INTERNAL;i++) lidlist[i]=nodes[i]->Lid();
333 return;
334 case BoundaryEnum:
335 for(int i=0;i<NUMNODES_BOUNDARY;i++) lidlist[i]=nodes[i]->Lid();
336 return;
337 default:
338 _error_("Numericalflux type " << EnumToStringx(type) << " not supported yet");
339 }
340}
341/*}}}*/
342/*FUNCTION Numericalflux::GetNumberOfNodes{{{*/
343int Numericalflux::GetNumberOfNodes(void){
344
345 int type;
346 inputs->GetInputValue(&type,NumericalfluxTypeEnum);
347
348 switch(type){
349 case InternalEnum:
350 return NUMNODES_INTERNAL;
351 case BoundaryEnum:
352 return NUMNODES_BOUNDARY;
353 default:
354 _error_("Numericalflux type " << EnumToStringx(type) << " not supported yet");
355 }
356
357}
358/*}}}*/
359/*FUNCTION Numericalflux::IsPenalty{{{*/
360bool Numericalflux::IsPenalty(void){
361 return false;
362}
363/*}}}*/
364/*FUNCTION Numericalflux::PenaltyCreateKMatrix {{{*/
365void Numericalflux::PenaltyCreateKMatrix(Matrix<IssmDouble>* Kff, Matrix<IssmDouble>* Kfs,IssmDouble kmax){
366
367 /*No stiffness loads applied, do nothing: */
368 return;
369
370}
371/*}}}*/
372/*FUNCTION Numericalflux::PenaltyCreatePVector{{{*/
373void Numericalflux::PenaltyCreatePVector(Vector<IssmDouble>* pf,IssmDouble kmax){
374
375 /*No penalty loads applied, do nothing: */
376 return;
377
378}
379/*}}}*/
380/*FUNCTION Numericalflux::InAnalysis{{{*/
381bool Numericalflux::InAnalysis(int in_analysis_type){
382 if (in_analysis_type==this->analysis_type) return true;
383 else return false;
384}
385/*}}}*/
386/*FUNCTION Numericalflux::SetwiseNodeConnectivity{{{*/
387void Numericalflux::SetwiseNodeConnectivity(int* pd_nz,int* po_nz,Node* node,bool* flags,int* flagsindices,int set1_enum,int set2_enum){
388
389 /*Output */
390 int d_nz = 0;
391 int o_nz = 0;
392
393 /*Loop over all nodes*/
394 for(int i=0;i<this->GetNumberOfNodes();i++){
395
396 if(!flags[this->nodes[i]->Lid()]){
397
398 /*flag current node so that no other element processes it*/
399 flags[this->nodes[i]->Lid()]=true;
400
401 int counter=0;
402 while(flagsindices[counter]>=0) counter++;
403 flagsindices[counter]=this->nodes[i]->Lid();
404
405 /*if node is clone, we have an off-diagonal non-zero, else it is a diagonal non-zero*/
406 switch(set2_enum){
407 case FsetEnum:
408 if(nodes[i]->indexing.fsize){
409 if(this->nodes[i]->IsClone())
410 o_nz += 1;
411 else
412 d_nz += 1;
413 }
414 break;
415 case GsetEnum:
416 if(nodes[i]->indexing.gsize){
417 if(this->nodes[i]->IsClone())
418 o_nz += 1;
419 else
420 d_nz += 1;
421 }
422 break;
423 case SsetEnum:
424 if(nodes[i]->indexing.ssize){
425 if(this->nodes[i]->IsClone())
426 o_nz += 1;
427 else
428 d_nz += 1;
429 }
430 break;
431 default: _error_("not supported");
432 }
433 }
434 }
435
436 /*Assign output pointers: */
437 *pd_nz=d_nz;
438 *po_nz=o_nz;
439}
440/*}}}*/
441
442/*Numericalflux management*/
443/*FUNCTION Numericalflux::CreateKMatrixMasstransport{{{*/
444ElementMatrix* Numericalflux::CreateKMatrixMasstransport(void){
445
446 int type;
447 inputs->GetInputValue(&type,NumericalfluxTypeEnum);
448
449 switch(type){
450 case InternalEnum:
451 return CreateKMatrixMasstransportInternal();
452 case BoundaryEnum:
453 return CreateKMatrixMasstransportBoundary();
454 default:
455 _error_("type not supported yet");
456 }
457}
458/*}}}*/
459/*FUNCTION Numericalflux::CreateKMatrixMasstransportInternal {{{*/
460ElementMatrix* Numericalflux::CreateKMatrixMasstransportInternal(void){
461
462 /* constants*/
463 const int numdof=NDOF1*NUMNODES_INTERNAL;
464
465 /* Intermediaries*/
466 int i,j,ig,index1,index2;
467 IssmDouble DL1,DL2,Jdet,dt,vx,vy,UdotN;
468 IssmDouble xyz_list[NUMVERTICES][3];
469 IssmDouble normal[2];
470 IssmDouble B[numdof];
471 IssmDouble Bprime[numdof];
472 IssmDouble Ke_g1[numdof][numdof];
473 IssmDouble Ke_g2[numdof][numdof];
474 GaussTria *gauss;
475
476 /*Initialize Element matrix and return if necessary*/
477 Tria* tria=(Tria*)element;
478 if(tria->NoIceInElement()) return NULL;
479 ElementMatrix* Ke=new ElementMatrix(nodes,NUMNODES_INTERNAL,this->parameters);
480
481 /*Retrieve all inputs and parameters*/
482 GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
483 parameters->FindParam(&dt,TimesteppingTimeStepEnum);
484 Input* vxaverage_input=tria->inputs->GetInput(VxEnum);
485 Input* vyaverage_input=tria->inputs->GetInput(VyEnum);
486 GetNormal(&normal[0],xyz_list);
487
488 /* Start looping on the number of gaussian points: */
489 index1=tria->GetNodeIndex(nodes[0]);
490 index2=tria->GetNodeIndex(nodes[1]);
491 gauss=new GaussTria(index1,index2,2);
492 for(ig=gauss->begin();ig<gauss->end();ig++){
493
494 gauss->GaussPoint(ig);
495
496 tria->GetSegmentBFlux(&B[0],gauss,index1,index2);
497 tria->GetSegmentBprimeFlux(&Bprime[0],gauss,index1,index2);
498
499 vxaverage_input->GetInputValue(&vx,gauss);
500 vyaverage_input->GetInputValue(&vy,gauss);
501 UdotN=vx*normal[0]+vy*normal[1];
502 tria->GetSegmentJacobianDeterminant(&Jdet,&xyz_list[0][0],gauss);
503 DL1=gauss->weight*Jdet*dt*UdotN/2;
504 DL2=gauss->weight*Jdet*dt*fabs(UdotN)/2;
505
506 TripleMultiply(&B[0],1,numdof,1,
507 &DL1,1,1,0,
508 &Bprime[0],1,numdof,0,
509 &Ke_g1[0][0],0);
510 TripleMultiply(&B[0],1,numdof,1,
511 &DL2,1,1,0,
512 &B[0],1,numdof,0,
513 &Ke_g2[0][0],0);
514
515 for(i=0;i<numdof;i++) for(j=0;j<numdof;j++) Ke->values[i*numdof+j]+=Ke_g1[i][j];
516 for(i=0;i<numdof;i++) for(j=0;j<numdof;j++) Ke->values[i*numdof+j]+=Ke_g2[i][j];
517 }
518
519 /*Clean up and return*/
520 delete gauss;
521 return Ke;
522}
523/*}}}*/
524/*FUNCTION Numericalflux::CreateKMatrixMasstransportBoundary {{{*/
525ElementMatrix* Numericalflux::CreateKMatrixMasstransportBoundary(void){
526
527 /* constants*/
528 const int numdof=NDOF1*NUMNODES_BOUNDARY;
529
530 /* Intermediaries*/
531 int i,j,ig,index1,index2;
532 IssmDouble DL,Jdet,dt,vx,vy,mean_vx,mean_vy,UdotN;
533 IssmDouble xyz_list[NUMVERTICES][3];
534 IssmDouble normal[2];
535 IssmDouble L[numdof];
536 IssmDouble Ke_g[numdof][numdof];
537 GaussTria *gauss;
538
539 /*Initialize Element matrix and return if necessary*/
540 ElementMatrix* Ke = NULL;
541 Tria* tria=(Tria*)element;
542 if(tria->NoIceInElement()) return NULL;
543
544 /*Retrieve all inputs and parameters*/
545 GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
546 parameters->FindParam(&dt,TimesteppingTimeStepEnum);
547 Input* vxaverage_input=tria->inputs->GetInput(VxEnum); _assert_(vxaverage_input);
548 Input* vyaverage_input=tria->inputs->GetInput(VyEnum); _assert_(vyaverage_input);
549 GetNormal(&normal[0],xyz_list);
550
551 /*Check wether it is an inflow or outflow BC (0 is the middle of the segment)*/
552 index1=tria->GetNodeIndex(nodes[0]);
553 index2=tria->GetNodeIndex(nodes[1]);
554
555 gauss=new GaussTria();
556 gauss->GaussEdgeCenter(index1,index2);
557 vxaverage_input->GetInputValue(&mean_vx,gauss);
558 vyaverage_input->GetInputValue(&mean_vy,gauss);
559 delete gauss;
560
561 UdotN=mean_vx*normal[0]+mean_vy*normal[1];
562 if (UdotN<=0){
563 return NULL; /*(u,n)<0 -> inflow, PenaltyCreatePVector will take care of it*/
564 }
565 else{
566 Ke=new ElementMatrix(nodes,NUMNODES_BOUNDARY,this->parameters);
567 }
568
569 /* Start looping on the number of gaussian points: */
570 gauss=new GaussTria(index1,index2,2);
571 for(ig=gauss->begin();ig<gauss->end();ig++){
572
573 gauss->GaussPoint(ig);
574
575 tria->GetSegmentNodalFunctions(&L[0],gauss,index1,index2);
576
577 vxaverage_input->GetInputValue(&vx,gauss);
578 vyaverage_input->GetInputValue(&vy,gauss);
579 UdotN=vx*normal[0]+vy*normal[1];
580 tria->GetSegmentJacobianDeterminant(&Jdet,&xyz_list[0][0],gauss);
581 DL=gauss->weight*Jdet*dt*UdotN;
582
583 TripleMultiply(&L[0],1,numdof,1,
584 &DL,1,1,0,
585 &L[0],1,numdof,0,
586 &Ke_g[0][0],0);
587
588 for(i=0;i<numdof;i++) for(j=0;j<numdof;j++) Ke->values[i*numdof+j]+=Ke_g[i][j];
589 }
590
591 /*Clean up and return*/
592 delete gauss;
593 return Ke;
594}
595/*}}}*/
596/*FUNCTION Numericalflux::CreateKMatrixBalancethickness{{{*/
597ElementMatrix* Numericalflux::CreateKMatrixBalancethickness(void){
598
599 int type;
600 inputs->GetInputValue(&type,NumericalfluxTypeEnum);
601
602 switch(type){
603 case InternalEnum:
604 return CreateKMatrixBalancethicknessInternal();
605 case BoundaryEnum:
606 return CreateKMatrixBalancethicknessBoundary();
607 default:
608 _error_("type not supported yet");
609 }
610}
611/*}}}*/
612/*FUNCTION Numericalflux::CreateKMatrixBalancethicknessInternal {{{*/
613ElementMatrix* Numericalflux::CreateKMatrixBalancethicknessInternal(void){
614
615 /* constants*/
616 const int numdof=NDOF1*NUMNODES_INTERNAL;
617
618 /* Intermediaries*/
619 int i,j,ig,index1,index2;
620 IssmDouble DL1,DL2,Jdet,vx,vy,UdotN;
621 IssmDouble xyz_list[NUMVERTICES][3];
622 IssmDouble normal[2];
623 IssmDouble B[numdof];
624 IssmDouble Bprime[numdof];
625 IssmDouble Ke_g1[numdof][numdof];
626 IssmDouble Ke_g2[numdof][numdof];
627 GaussTria *gauss;
628
629 /*Initialize Element matrix and return if necessary*/
630 Tria* tria=(Tria*)element;
631 if(tria->NoIceInElement()) return NULL;
632 ElementMatrix* Ke=new ElementMatrix(nodes,NUMNODES_INTERNAL,this->parameters);
633
634 /*Retrieve all inputs and parameters*/
635 GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
636 Input* vxaverage_input=tria->inputs->GetInput(VxEnum);
637 Input* vyaverage_input=tria->inputs->GetInput(VyEnum);
638 GetNormal(&normal[0],xyz_list);
639
640 /* Start looping on the number of gaussian points: */
641 index1=tria->GetNodeIndex(nodes[0]);
642 index2=tria->GetNodeIndex(nodes[1]);
643 gauss=new GaussTria(index1,index2,2);
644 for(ig=gauss->begin();ig<gauss->end();ig++){
645
646 gauss->GaussPoint(ig);
647
648 tria->GetSegmentBFlux(&B[0],gauss,index1,index2);
649 tria->GetSegmentBprimeFlux(&Bprime[0],gauss,index1,index2);
650
651 vxaverage_input->GetInputValue(&vx,gauss);
652 vyaverage_input->GetInputValue(&vy,gauss);
653 UdotN=vx*normal[0]+vy*normal[1];
654 tria->GetSegmentJacobianDeterminant(&Jdet,&xyz_list[0][0],gauss);
655 DL1=gauss->weight*Jdet*UdotN/2;
656 DL2=gauss->weight*Jdet*fabs(UdotN)/2;
657
658 TripleMultiply(&B[0],1,numdof,1,
659 &DL1,1,1,0,
660 &Bprime[0],1,numdof,0,
661 &Ke_g1[0][0],0);
662 TripleMultiply(&B[0],1,numdof,1,
663 &DL2,1,1,0,
664 &B[0],1,numdof,0,
665 &Ke_g2[0][0],0);
666
667 for(i=0;i<numdof;i++) for(j=0;j<numdof;j++) Ke->values[i*numdof+j]+=Ke_g1[i][j];
668 for(i=0;i<numdof;i++) for(j=0;j<numdof;j++) Ke->values[i*numdof+j]+=Ke_g2[i][j];
669 }
670
671 /*Clean up and return*/
672 delete gauss;
673 return Ke;
674}
675/*}}}*/
676/*FUNCTION Numericalflux::CreateKMatrixBalancethicknessBoundary {{{*/
677ElementMatrix* Numericalflux::CreateKMatrixBalancethicknessBoundary(void){
678
679 /* constants*/
680 const int numdof=NDOF1*NUMNODES_BOUNDARY;
681
682 /* Intermediaries*/
683 int i,j,ig,index1,index2;
684 IssmDouble DL,Jdet,vx,vy,mean_vx,mean_vy,UdotN;
685 IssmDouble xyz_list[NUMVERTICES][3];
686 IssmDouble normal[2];
687 IssmDouble L[numdof];
688 IssmDouble Ke_g[numdof][numdof];
689 GaussTria *gauss;
690
691 /*Initialize Element matrix and return if necessary*/
692 ElementMatrix* Ke = NULL;
693 Tria* tria=(Tria*)element;
694 if(tria->NoIceInElement()) return NULL;
695
696 /*Retrieve all inputs and parameters*/
697 GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
698 Input* vxaverage_input=tria->inputs->GetInput(VxEnum);
699 Input* vyaverage_input=tria->inputs->GetInput(VyEnum);
700 GetNormal(&normal[0],xyz_list);
701
702 /*Check wether it is an inflow or outflow BC (0 is the middle of the segment)*/
703 index1=tria->GetNodeIndex(nodes[0]);
704 index2=tria->GetNodeIndex(nodes[1]);
705
706 gauss=new GaussTria();
707 gauss->GaussEdgeCenter(index1,index2);
708 vxaverage_input->GetInputValue(&mean_vx,gauss);
709 vyaverage_input->GetInputValue(&mean_vy,gauss);
710 delete gauss;
711
712 UdotN=mean_vx*normal[0]+mean_vy*normal[1];
713 if (UdotN<=0){
714 return NULL; /*(u,n)<0 -> inflow, PenaltyCreatePVector will take care of it*/
715 }
716 else{
717 Ke=new ElementMatrix(nodes,NUMNODES_BOUNDARY,this->parameters);
718 }
719
720 /* Start looping on the number of gaussian points: */
721 gauss=new GaussTria(index1,index2,2);
722 for(ig=gauss->begin();ig<gauss->end();ig++){
723
724 gauss->GaussPoint(ig);
725
726 tria->GetSegmentNodalFunctions(&L[0],gauss,index1,index2);
727
728 vxaverage_input->GetInputValue(&vx,gauss);
729 vyaverage_input->GetInputValue(&vy,gauss);
730 UdotN=vx*normal[0]+vy*normal[1];
731 tria->GetSegmentJacobianDeterminant(&Jdet,&xyz_list[0][0],gauss);
732 DL=gauss->weight*Jdet*UdotN;
733
734 TripleMultiply(&L[0],1,numdof,1,
735 &DL,1,1,0,
736 &L[0],1,numdof,0,
737 &Ke_g[0][0],0);
738
739 for(i=0;i<numdof;i++) for(j=0;j<numdof;j++) Ke->values[i*numdof+j]+=Ke_g[i][j];
740 }
741
742 /*Clean up and return*/
743 delete gauss;
744 return Ke;
745}
746/*}}}*/
747/*FUNCTION Numericalflux::CreateKMatrixAdjointBalancethickness{{{*/
748ElementMatrix* Numericalflux::CreateKMatrixAdjointBalancethickness(void){
749
750 int type;
751 inputs->GetInputValue(&type,NumericalfluxTypeEnum);
752
753 switch(type){
754 case InternalEnum:
755 return CreateKMatrixAdjointBalancethicknessInternal();
756 case BoundaryEnum:
757 return CreateKMatrixAdjointBalancethicknessBoundary();
758 default:
759 _error_("type not supported yet");
760 }
761}
762/*}}}*/
763/*FUNCTION Numericalflux::CreateKMatrixAdjointBalancethicknessInternal {{{*/
764ElementMatrix* Numericalflux::CreateKMatrixAdjointBalancethicknessInternal(void){
765
766 ElementMatrix* Ke=CreateKMatrixBalancethicknessInternal();
767 if (Ke) Ke->Transpose();
768 return Ke;
769}
770/*}}}*/
771/*FUNCTION Numericalflux::CreateKMatrixAdjointBalancethicknessBoundary {{{*/
772ElementMatrix* Numericalflux::CreateKMatrixAdjointBalancethicknessBoundary(void){
773
774 ElementMatrix* Ke=CreateKMatrixBalancethicknessBoundary();
775 if(Ke) Ke->Transpose();
776 return Ke;
777}
778/*}}}*/
779/*FUNCTION Numericalflux::CreatePVectorMasstransport{{{*/
780ElementVector* Numericalflux::CreatePVectorMasstransport(void){
781
782 int type;
783 inputs->GetInputValue(&type,NumericalfluxTypeEnum);
784
785 switch(type){
786 case InternalEnum:
787 return CreatePVectorMasstransportInternal();
788 case BoundaryEnum:
789 return CreatePVectorMasstransportBoundary();
790 default:
791 _error_("type not supported yet");
792 }
793}
794/*}}}*/
795/*FUNCTION Numericalflux::CreatePVectorMasstransportInternal{{{*/
796ElementVector* Numericalflux::CreatePVectorMasstransportInternal(void){
797
798 /*Nothing added to PVector*/
799 return NULL;
800
801}
802/*}}}*/
803/*FUNCTION Numericalflux::CreatePVectorMasstransportBoundary{{{*/
804ElementVector* Numericalflux::CreatePVectorMasstransportBoundary(void){
805
806 /* constants*/
807 const int numdof=NDOF1*NUMNODES_BOUNDARY;
808
809 /* Intermediaries*/
810 int i,ig,index1,index2;
811 IssmDouble DL,Jdet,dt,vx,vy,mean_vx,mean_vy,UdotN,thickness;
812 IssmDouble xyz_list[NUMVERTICES][3];
813 IssmDouble normal[2];
814 IssmDouble L[numdof];
815 GaussTria *gauss;
816
817 /*Initialize Load Vector and return if necessary*/
818 ElementVector* pe = NULL;
819 Tria* tria=(Tria*)element;
820 if(tria->NoIceInElement()) return NULL;
821
822 /*Retrieve all inputs and parameters*/
823 GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
824 parameters->FindParam(&dt,TimesteppingTimeStepEnum);
825 Input* vxaverage_input =tria->inputs->GetInput(VxEnum); _assert_(vxaverage_input);
826 Input* vyaverage_input =tria->inputs->GetInput(VyEnum); _assert_(vyaverage_input);
827 Input* spcthickness_input=tria->inputs->GetInput(MasstransportSpcthicknessEnum); _assert_(spcthickness_input);
828 GetNormal(&normal[0],xyz_list);
829
830 /*Check wether it is an inflow or outflow BC (0 is the middle of the segment)*/
831 index1=tria->GetNodeIndex(nodes[0]);
832 index2=tria->GetNodeIndex(nodes[1]);
833
834 gauss=new GaussTria();
835 gauss->GaussEdgeCenter(index1,index2);
836 vxaverage_input->GetInputValue(&mean_vx,gauss);
837 vyaverage_input->GetInputValue(&mean_vy,gauss);
838 delete gauss;
839
840 UdotN=mean_vx*normal[0]+mean_vy*normal[1];
841 if (UdotN>0){
842 return NULL; /*(u,n)>0 -> outflow, PenaltyCreateKMatrix will take care of it*/
843 }
844 else{
845 pe=new ElementVector(nodes,NUMNODES_BOUNDARY,this->parameters);
846 }
847
848 /* Start looping on the number of gaussian points: */
849 gauss=new GaussTria(index1,index2,2);
850 for(ig=gauss->begin();ig<gauss->end();ig++){
851
852 gauss->GaussPoint(ig);
853
854 tria->GetSegmentNodalFunctions(&L[0],gauss,index1,index2);
855
856 vxaverage_input->GetInputValue(&vx,gauss);
857 vyaverage_input->GetInputValue(&vy,gauss);
858 spcthickness_input->GetInputValue(&thickness,gauss);
859 if(xIsNan<IssmDouble>(thickness)) _error_("Cannot weakly apply constraint because NaN was provided");
860
861 UdotN=vx*normal[0]+vy*normal[1];
862 tria->GetSegmentJacobianDeterminant(&Jdet,&xyz_list[0][0],gauss);
863 DL= - gauss->weight*Jdet*dt*UdotN*thickness;
864
865 for(i=0;i<numdof;i++) pe->values[i] += DL*L[i];
866 }
867
868 /*Clean up and return*/
869 delete gauss;
870 return pe;
871}
872/*}}}*/
873/*FUNCTION Numericalflux::CreatePVectorBalancethickness{{{*/
874ElementVector* Numericalflux::CreatePVectorBalancethickness(void){
875
876 int type;
877 inputs->GetInputValue(&type,NumericalfluxTypeEnum);
878
879 switch(type){
880 case InternalEnum:
881 return CreatePVectorBalancethicknessInternal();
882 case BoundaryEnum:
883 return CreatePVectorBalancethicknessBoundary();
884 default:
885 _error_("type not supported yet");
886 }
887}
888/*}}}*/
889/*FUNCTION Numericalflux::CreatePVectorBalancethicknessInternal{{{*/
890ElementVector* Numericalflux::CreatePVectorBalancethicknessInternal(void){
891
892 /*Nothing added to PVector*/
893 return NULL;
894
895}
896/*}}}*/
897/*FUNCTION Numericalflux::CreatePVectorBalancethicknessBoundary{{{*/
898ElementVector* Numericalflux::CreatePVectorBalancethicknessBoundary(void){
899
900 /* constants*/
901 const int numdof=NDOF1*NUMNODES_BOUNDARY;
902
903 /* Intermediaries*/
904 int i,ig,index1,index2;
905 IssmDouble DL,Jdet,vx,vy,mean_vx,mean_vy,UdotN,thickness;
906 IssmDouble xyz_list[NUMVERTICES][3];
907 IssmDouble normal[2];
908 IssmDouble L[numdof];
909 GaussTria *gauss;
910
911 /*Initialize Load Vector and return if necessary*/
912 ElementVector* pe = NULL;
913 Tria* tria=(Tria*)element;
914 if(tria->NoIceInElement()) return NULL;
915
916 /*Retrieve all inputs and parameters*/
917 GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
918 Input* vxaverage_input=tria->inputs->GetInput(VxEnum); _assert_(vxaverage_input);
919 Input* vyaverage_input=tria->inputs->GetInput(VyEnum); _assert_(vyaverage_input);
920 Input* thickness_input=tria->inputs->GetInput(ThicknessEnum); _assert_(thickness_input);
921 GetNormal(&normal[0],xyz_list);
922
923 /*Check wether it is an inflow or outflow BC (0 is the middle of the segment)*/
924 index1=tria->GetNodeIndex(nodes[0]);
925 index2=tria->GetNodeIndex(nodes[1]);
926
927 gauss=new GaussTria();
928 gauss->GaussEdgeCenter(index1,index2);
929 vxaverage_input->GetInputValue(&mean_vx,gauss);
930 vyaverage_input->GetInputValue(&mean_vy,gauss);
931 delete gauss;
932 UdotN=mean_vx*normal[0]+mean_vy*normal[1];
933 if (UdotN>0){
934 return NULL; /*(u,n)>0 -> outflow, PenaltyCreateKMatrix will take care of it*/
935 }
936 else{
937 pe=new ElementVector(nodes,NUMNODES_BOUNDARY,this->parameters);
938 }
939
940 /* Start looping on the number of gaussian points: */
941 gauss=new GaussTria(index1,index2,2);
942 for(ig=gauss->begin();ig<gauss->end();ig++){
943
944 gauss->GaussPoint(ig);
945
946 tria->GetSegmentNodalFunctions(&L[0],gauss,index1,index2);
947
948 vxaverage_input->GetInputValue(&vx,gauss);
949 vyaverage_input->GetInputValue(&vy,gauss);
950 thickness_input->GetInputValue(&thickness,gauss);
951
952 UdotN=vx*normal[0]+vy*normal[1];
953 tria->GetSegmentJacobianDeterminant(&Jdet,&xyz_list[0][0],gauss);
954 DL= - gauss->weight*Jdet*UdotN*thickness;
955
956 for(i=0;i<numdof;i++) pe->values[i] += DL*L[i];
957 }
958
959 /*Clean up and return*/
960 delete gauss;
961 return pe;
962}
963/*}}}*/
964/*FUNCTION Numericalflux::CreatePVectorAdjointBalancethickness{{{*/
965ElementVector* Numericalflux::CreatePVectorAdjointBalancethickness(void){
966
967 /*No PVector for the Adjoint*/
968 return NULL;
969}
970/*}}}*/
971/*FUNCTION Numericalflux::GetNormal {{{*/
972void Numericalflux:: GetNormal(IssmDouble* normal,IssmDouble xyz_list[4][3]){
973
974 /*Build unit outward pointing vector*/
975 IssmDouble vector[2];
976 IssmDouble norm;
977
978 vector[0]=xyz_list[1][0] - xyz_list[0][0];
979 vector[1]=xyz_list[1][1] - xyz_list[0][1];
980
981 norm=sqrt(pow(vector[0],2.0)+pow(vector[1],2.0));
982
983 normal[0]= + vector[1]/norm;
984 normal[1]= - vector[0]/norm;
985}
986/*}}}*/
Note: See TracBrowser for help on using the repository browser.