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

Last change on this file since 15396 was 15373, checked in by Mathieu Morlighem, 12 years ago

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