source: issm/branches/trunk-jpl-damage/src/c/classes/objects/Loads/Numericalflux.cpp@ 12878

Last change on this file since 12878 was 12832, checked in by Eric.Larour, 13 years ago

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