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

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

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