source: issm/branches/trunk-jpl-damage/src/c/objects/Node.cpp@ 11427

Last change on this file since 11427 was 11001, checked in by Mathieu Morlighem, 13 years ago

Results on vertices are now averaged using node connectivity as in PatchToVec
File size: 28.4 KB
Line 
1/*!\file Node.c
2 * \brief: implementation of the Node object
3 */
4
5/*Include files: {{{1*/
6#ifdef HAVE_CONFIG_H
7 #include <config.h>
8#else
9#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!"
10#endif
11
12#include <stdio.h>
13#include <string.h>
14#include "./objects.h"
15#include "../Container/Container.h"
16#include "../EnumDefinitions/EnumDefinitions.h"
17#include "../shared/shared.h"
18#include "../include/include.h"
19#include "../modules/modules.h"
20/*}}}*/
21
22/*Node constructors and destructors:*/
23/*FUNCTION Node::Node() default constructor {{{1*/
24Node::Node(){
25 this->inputs=NULL;
26 this->hvertex=NULL;
27 return;
28}
29/*}}}*/
30/*FUNCTION Node::Node(int node_id,int node_sid,int vertex_id,int io_index, IoModel* iomodel,int analysis_type) {{{1*/
31Node::Node(int node_id,int node_sid,int vertex_id,int io_index, IoModel* iomodel,int analysis_type){
32
33 /*Intermediary*/
34 int k,l;
35 int gsize;
36 int dim;
37
38 /*Fetch parameters: */
39 iomodel->Constant(&dim,MeshDimensionEnum);
40
41 /*id: */
42 this->id=node_id;
43 this->sid=node_sid;
44 this->analysis_type=analysis_type;
45
46 /*Initialize coord_system: Identity matrix by default*/
47 for(k=0;k<3;k++) for(l=0;l<3;l++) this->coord_system[k][l]=0.0;
48 for(k=0;k<3;k++) this->coord_system[k][k]=1.0;
49
50 /*indexing:*/
51 DistributeNumDofs(&this->indexing,analysis_type,iomodel->Data(FlowequationVertexEquationEnum)+io_index); //number of dofs per node
52 gsize=this->indexing.gsize;
53
54 /*Hooks*/
55 this->hvertex=new Hook(&vertex_id,1); //node id is the same as the vertex id, continuous galerkin!
56
57 //intialize inputs, and add as many inputs per element as requested:
58 this->inputs=new Inputs();
59 if (iomodel->Data(MeshVertexonbedEnum))
60 this->inputs->AddInput(new BoolInput(MeshVertexonbedEnum,(IssmBool)iomodel->Data(MeshVertexonbedEnum)[io_index]));
61 if (iomodel->Data(MeshVertexonsurfaceEnum))
62 this->inputs->AddInput(new BoolInput(MeshVertexonsurfaceEnum,(IssmBool)iomodel->Data(MeshVertexonsurfaceEnum)[io_index]));
63 if (iomodel->Data(MaskVertexonfloatingiceEnum))
64 this->inputs->AddInput(new BoolInput(MaskVertexonfloatingiceEnum,(IssmBool)iomodel->Data(MaskVertexonfloatingiceEnum)[io_index]));
65 if (iomodel->Data(MaskVertexongroundediceEnum))
66 this->inputs->AddInput(new BoolInput(MaskVertexongroundediceEnum,(IssmBool)iomodel->Data(MaskVertexongroundediceEnum)[io_index]));
67 if (analysis_type==DiagnosticHorizAnalysisEnum)
68 this->inputs->AddInput(new IntInput(ApproximationEnum,(IssmInt)iomodel->Data(FlowequationVertexEquationEnum)[io_index]));
69
70 /*set single point constraints: */
71
72 /*spc all nodes on water*/
73 if (!iomodel->Data(MaskVertexonwaterEnum)) _error_("iomodel->nodeonwater is NULL");
74 if (iomodel->Data(MaskVertexonwaterEnum)[io_index]){
75 for(k=1;k<=gsize;k++){
76 this->FreezeDof(k);
77 }
78 }
79
80 /*Diagnostic Horiz*/
81 #ifdef _HAVE_DIAGNOSTIC_
82 if (analysis_type==DiagnosticHorizAnalysisEnum){
83
84 /*Coordinate system provided, convert to coord_system matrix*/
85 _assert_(iomodel->Data(DiagnosticReferentialEnum));
86 XZvectorsToCoordinateSystem(&this->coord_system[0][0],iomodel->Data(DiagnosticReferentialEnum)+io_index*6);
87
88 if (dim==3){
89 /*We have a 3d mesh, we may have collapsed elements, hence dead nodes. Freeze them out: */
90 _assert_(iomodel->Data(MeshVertexonbedEnum));
91 _assert_(iomodel->Data(FlowequationVertexEquationEnum));
92 if (iomodel->Data(FlowequationVertexEquationEnum)[io_index]==MacAyealApproximationEnum && !iomodel->Data(MeshVertexonbedEnum)[io_index]){
93 for(k=1;k<=gsize;k++) this->FreezeDof(k);
94 }
95 if (iomodel->Data(FlowequationVertexEquationEnum)[io_index]==MacAyealPattynApproximationEnum && iomodel->Data(FlowequationBordermacayealEnum)[io_index]){
96 if(!iomodel->Data(MeshVertexonbedEnum)[io_index]){
97 for(k=1;k<=gsize;k++) this->FreezeDof(k);
98 }
99 }
100 if (iomodel->Data(FlowequationVertexEquationEnum)[io_index]==MacAyealStokesApproximationEnum && iomodel->Data(FlowequationBordermacayealEnum)[io_index]){
101 if(!iomodel->Data(MeshVertexonbedEnum)[io_index]){
102 for(k=1;k<=2;k++) this->FreezeDof(k);
103 }
104 }
105 }
106 /*spc all nodes on hutter*/
107 if (iomodel->Data(FlowequationVertexEquationEnum)[io_index]==HutterApproximationEnum){
108 for(k=1;k<=gsize;k++){
109 this->FreezeDof(k);
110 }
111 }
112 }
113 #endif
114
115 /*Diagnostic Hutter*/
116 if (analysis_type==DiagnosticHutterAnalysisEnum){
117 if (!iomodel->Data(FlowequationVertexEquationEnum)) _error_("iomodel->vertices_type is NULL");
118 /*Constrain all nodes that are not Hutter*/
119 if (!iomodel->Data(FlowequationVertexEquationEnum)[io_index]==HutterApproximationEnum){
120 for(k=1;k<=gsize;k++){
121 this->FreezeDof(k);
122 }
123 }
124 }
125
126 /*Prognostic/ Melting/ Slopecompute/ Balancethickness*/
127 if (
128 analysis_type==PrognosticAnalysisEnum ||
129 analysis_type==MeltingAnalysisEnum ||
130 analysis_type==BedSlopeAnalysisEnum ||
131 analysis_type==SurfaceSlopeAnalysisEnum ||
132 analysis_type==BalancethicknessAnalysisEnum
133 ){
134 if (dim==3){
135 /*On a 3d mesh, we may have collapsed elements, hence dead nodes. Freeze them out: */
136 _assert_(iomodel->Data(MeshVertexonbedEnum));
137 if (!iomodel->Data(MeshVertexonbedEnum)[io_index]){
138 for(k=1;k<=gsize;k++){
139 this->FreezeDof(k);
140 }
141 }
142 }
143 }
144
145}
146/*}}}*/
147/*FUNCTION Node::~Node(){{{1*/
148Node::~Node(){
149 delete inputs;
150 delete hvertex;
151 return;
152}
153/*}}}*/
154
155/*Object virtual functions definitions:*/
156/*FUNCTION Node::Echo{{{1*/
157void Node::Echo(void){
158
159 printf("Node:\n");
160 printf(" id: %i\n",id);
161 printf(" sid: %i\n",sid);
162 printf(" analysis_type: %s\n",EnumToStringx(analysis_type));
163 indexing.Echo();
164 printf(" hvertex: not displayed\n");
165 printf(" inputs: %p\n",inputs);
166
167
168}
169/*}}}*/
170/*FUNCTION Node::DeepEcho{{{1*/
171void Node::DeepEcho(void){
172
173 printf("Node:\n");
174 printf(" id: %i\n",id);
175 printf(" sid: %i\n",sid);
176 printf(" analysis_type: %s\n",EnumToStringx(analysis_type));
177 indexing.DeepEcho();
178 printf("Vertex:\n");
179 hvertex->DeepEcho();
180 printf(" inputs\n");
181
182
183}
184/*}}}*/
185/*FUNCTION Node::Id{{{1*/
186int Node::Id(void){ return id; }
187/*}}}*/
188/*FUNCTION Node::MyRank{{{1*/
189int Node::MyRank(void){
190 extern int my_rank;
191
192 return my_rank;
193}
194/*}}}*/
195#ifdef _SERIAL_
196/*FUNCTION Node::Marshall{{{1*/
197void Node::Marshall(char** pmarshalled_dataset){
198
199 char* marshalled_dataset=NULL;
200 int enum_type=0;
201 char* marshalled_inputs=NULL;
202 int marshalled_inputssize;
203
204 /*recover marshalled_dataset: */
205 marshalled_dataset=*pmarshalled_dataset;
206
207 /*get enum type of Node: */
208 enum_type=NodeEnum;
209
210 /*marshall enum: */
211 memcpy(marshalled_dataset,&enum_type,sizeof(enum_type));marshalled_dataset+=sizeof(enum_type);
212
213 /*marshall Node data: */
214 memcpy(marshalled_dataset,&id,sizeof(id));marshalled_dataset+=sizeof(id);
215 memcpy(marshalled_dataset,&sid,sizeof(sid));marshalled_dataset+=sizeof(sid);
216 memcpy(marshalled_dataset,&analysis_type,sizeof(analysis_type));marshalled_dataset+=sizeof(analysis_type);
217 memcpy(marshalled_dataset,&coord_system,9*sizeof(double));marshalled_dataset+=9*sizeof(double);
218
219 /*marshall objects: */
220 indexing.Marshall(&marshalled_dataset);
221 hvertex->Marshall(&marshalled_dataset);
222
223 /*Marshall inputs: */
224 marshalled_inputssize=inputs->MarshallSize();
225 marshalled_inputs=inputs->Marshall();
226 memcpy(marshalled_dataset,marshalled_inputs,marshalled_inputssize*sizeof(char));
227 marshalled_dataset+=marshalled_inputssize;
228
229 /*Free ressources:*/
230 xfree((void**)&marshalled_inputs);
231
232 *pmarshalled_dataset=marshalled_dataset;
233 return;
234}
235/*}}}*/
236/*FUNCTION Node::MarshallSize{{{1*/
237int Node::MarshallSize(){
238
239 return sizeof(id)+
240 sizeof(sid)+
241 indexing.MarshallSize()+
242 hvertex->MarshallSize()+
243 inputs->MarshallSize()+
244 sizeof(analysis_type)+
245 9*sizeof(double)+
246 sizeof(int); //sizeof(int) for enum type
247}
248/*}}}*/
249/*FUNCTION Node::Demarshall{{{1*/
250void Node::Demarshall(char** pmarshalled_dataset){
251
252 char* marshalled_dataset=NULL;
253
254 /*recover marshalled_dataset: */
255 marshalled_dataset=*pmarshalled_dataset;
256
257 /*this time, no need to get enum type, the pointer directly points to the beginning of the
258 *object data (thanks to DataSet::Demarshall):*/
259
260 memcpy(&id,marshalled_dataset,sizeof(id));marshalled_dataset+=sizeof(id);
261 memcpy(&sid,marshalled_dataset,sizeof(sid));marshalled_dataset+=sizeof(sid);
262 memcpy(&analysis_type,marshalled_dataset,sizeof(analysis_type));marshalled_dataset+=sizeof(analysis_type);
263 memcpy(&coord_system,marshalled_dataset,9*sizeof(double));marshalled_dataset+=9*sizeof(double);
264
265 /*demarshall objects: */
266 indexing.Demarshall(&marshalled_dataset);
267 hvertex=new Hook(); hvertex->Demarshall(&marshalled_dataset);
268
269 /*demarshall inputs: */
270 inputs=(Inputs*)DataSetDemarshallRaw(&marshalled_dataset);
271
272 /*return: */
273 *pmarshalled_dataset=marshalled_dataset;
274 return;
275}
276/*}}}*/
277#endif
278/*FUNCTION Node::ObjectEnum{{{1*/
279int Node::ObjectEnum(void){
280
281 return NodeEnum;
282
283}
284/*}}}*/
285
286/*Node management:*/
287/*FUNCTION Node::Configure {{{1*/
288void Node::Configure(DataSet* nodesin,Vertices* verticesin){
289
290 int i;
291
292 /*Take care of hooking up all objects for this element, ie links the objects in the hooks to their respective
293 * datasets, using internal ids and off_sets hidden in hooks: */
294 hvertex->configure(verticesin);
295
296}
297/*FUNCTION Node::SetCurrentConfiguration {{{1*/
298void Node::SetCurrentConfiguration(DataSet* nodesin,Vertices* verticesin){
299
300}
301/*FUNCTION Node::GetDof {{{1*/
302int Node::GetDof(int dofindex,int setenum){
303
304 if(setenum==GsetEnum){
305 _assert_(dofindex>=0 && dofindex<indexing.gsize);
306 return indexing.gdoflist[dofindex];
307 }
308 else if(setenum==FsetEnum){
309 _assert_(dofindex>=0 && dofindex<indexing.fsize);
310 return indexing.fdoflist[dofindex];
311 }
312 else if(setenum==SsetEnum){
313 _assert_(dofindex>=0 && dofindex<indexing.ssize);
314 return indexing.sdoflist[dofindex];
315 }
316 else _error_("%s%s%s"," set of enum type ",EnumToStringx(setenum)," not supported yet!");
317
318}
319/*}}}*/
320/*FUNCTION Node::GetDofList1{{{1*/
321int Node::GetDofList1(void){
322
323 Vertex* vertex=NULL;
324
325 vertex=(Vertex*)this->hvertex->delivers();
326
327 return vertex->dof;
328}
329/*}}}*/
330/*FUNCTION Node::GetDofList{{{1*/
331void Node::GetDofList(int* outdoflist,int approximation_enum,int setenum){
332 int i;
333 int count=0;
334 int count2=0;
335
336 if(approximation_enum==NoneApproximationEnum){
337 if(setenum==GsetEnum)for(i=0;i<this->indexing.gsize;i++) outdoflist[i]=indexing.gdoflist[i];
338 if(setenum==FsetEnum)for(i=0;i<this->indexing.fsize;i++) outdoflist[i]=indexing.fdoflist[i];
339 if(setenum==SsetEnum)for(i=0;i<this->indexing.ssize;i++) outdoflist[i]=indexing.sdoflist[i];
340 }
341 else{
342
343 if(setenum==GsetEnum){
344 if(indexing.doftype){
345 count=0;
346 for(i=0;i<this->indexing.gsize;i++){
347 if(indexing.doftype[i]==approximation_enum){
348 outdoflist[count]=indexing.gdoflist[i];
349 count++;
350 }
351 }
352 _assert_(count); //at least one dof should be the approximation requested
353 }
354 else for(i=0;i<this->indexing.gsize;i++) outdoflist[i]=indexing.gdoflist[i];
355 }
356 else if(setenum==FsetEnum){
357 if(indexing.doftype){
358 count=0;
359 count2=0;
360 for(i=0;i<this->indexing.gsize;i++){
361 if(indexing.f_set[i]){
362 if(indexing.doftype[i]==approximation_enum){
363 outdoflist[count]=indexing.fdoflist[count2];
364 count++;
365 }
366 count2++;
367 }
368 }
369 }
370 else for(i=0;i<this->indexing.fsize;i++) outdoflist[i]=indexing.fdoflist[i];
371 }
372 else if(setenum==SsetEnum){
373 if(indexing.doftype){
374 count=0;
375 count2=0;
376 for(i=0;i<this->indexing.gsize;i++){
377 if(indexing.s_set[i]){
378 if(indexing.doftype[i]==approximation_enum){
379 outdoflist[count]=indexing.sdoflist[count2];
380 count++;
381 }
382 count2++;
383 }
384 }
385 }
386 else for(i=0;i<this->indexing.ssize;i++) outdoflist[i]=indexing.sdoflist[i];
387 }
388 else _error_("%s%s%s"," set of enum type ",EnumToStringx(setenum)," not supported yet!");
389 }
390}
391/*}}}*/
392/*FUNCTION Node::GetSidList{{{1*/
393int Node::GetSidList(void){
394
395 Vertex* vertex=NULL;
396
397 vertex=(Vertex*)this->hvertex->delivers();
398
399 return vertex->sid;
400}
401/*}}}*/
402/*FUNCTION Node::GetLocalDofList{{{1*/
403void Node::GetLocalDofList(int* outdoflist,int approximation_enum,int setenum){
404 int i;
405 int count=0;
406 int count2=0;
407
408 if(approximation_enum==NoneApproximationEnum){
409 if(setenum==GsetEnum)for(i=0;i<this->indexing.gsize;i++) outdoflist[i]=i;
410 else if(setenum==FsetEnum){
411 count=0;
412 for(i=0;i<this->indexing.gsize;i++){
413 if(indexing.f_set[i]){
414 outdoflist[count]=i;
415 count++;
416 }
417 }
418 }
419 else if(setenum==SsetEnum){
420 count=0;
421 for(i=0;i<this->indexing.gsize;i++){
422 if(indexing.s_set[i]){
423 outdoflist[count]=i;
424 count++;
425 }
426 }
427 }
428 else _error_("%s%s%s"," set of enum type ",EnumToStringx(setenum)," not supported yet!");
429 }
430 else{
431
432 if(setenum==GsetEnum){
433 if(indexing.doftype){
434 count=0;
435 for(i=0;i<this->indexing.gsize;i++){
436 if(indexing.doftype[i]==approximation_enum){
437 outdoflist[count]=count;
438 count++;
439 }
440 }
441 _assert_(count);
442 }
443 else for(i=0;i<this->indexing.gsize;i++) outdoflist[i]=i;
444 }
445 else if(setenum==FsetEnum){
446
447 if(indexing.doftype){
448 count=0;
449 count2=0;
450 for(i=0;i<this->indexing.gsize;i++){
451 if(indexing.doftype[i]==approximation_enum){
452 if(indexing.f_set[i]){
453 outdoflist[count]=count2;
454 count++;
455 }
456 count2++;
457 }
458 }
459 _assert_(count2);
460 }
461 else{
462
463 count=0;
464 for(i=0;i<this->indexing.gsize;i++){
465 if(indexing.f_set[i]){
466 outdoflist[count]=i;
467 count++;
468 }
469 }
470 }
471 }
472 else if(setenum==SsetEnum){
473 if(indexing.doftype){
474 count=0;
475 count2=0;
476 for(i=0;i<this->indexing.gsize;i++){
477 if(indexing.doftype[i]==approximation_enum){
478 if(indexing.s_set[i]){
479 outdoflist[count]=count2;
480 count++;
481 }
482 count2++;
483 }
484 }
485 _assert_(count2);
486 }
487 else{
488 count=0;
489 for(i=0;i<this->indexing.gsize;i++){
490 if(indexing.s_set[i]){
491 outdoflist[count]=i;
492 count++;
493 }
494 }
495 }
496 }
497 else _error_("%s%s%s"," set of enum type ",EnumToStringx(setenum)," not supported yet!");
498 }
499}
500/*}}}*/
501/*FUNCTION Node::Sid{{{1*/
502int Node::Sid(void){ return sid; }
503/*}}}*/
504/*FUNCTION Node::GetVertexId {{{1*/
505int Node::GetVertexId(void){
506
507 Vertex* vertex=NULL;
508
509 vertex=(Vertex*)hvertex->delivers();
510 return vertex->id;
511}
512/*}}}*/
513/*FUNCTION Node::GetVertexDof {{{1*/
514int Node::GetVertexDof(void){
515
516 Vertex* vertex=NULL;
517
518 vertex=(Vertex*)hvertex->delivers();
519 return vertex->dof;
520}
521/*}}}*/
522#ifdef _HAVE_DIAGNOSTIC_
523/*FUNCTION Node::GetCoordinateSystem{{{1*/
524void Node::GetCoordinateSystem(double* coord_system_out){
525
526 /*Copy coord_system*/
527 for(int k=0;k<3;k++) for(int l=0;l<3;l++) coord_system_out[3*k+l]=this->coord_system[k][l];
528
529}
530/*}}}*/
531#endif
532/*FUNCTION Node::SetVertexDof {{{1*/
533void Node::SetVertexDof(int in_dof){
534
535 Vertex* vertex=NULL;
536
537 vertex=(Vertex*)hvertex->delivers();
538 vertex->dof=in_dof;
539
540}
541/*}}}*/
542/*FUNCTION Node::InAnalysis{{{1*/
543bool Node::InAnalysis(int in_analysis_type){
544 if (in_analysis_type==this->analysis_type) return true;
545 else return false;
546}
547/*}}}*/
548
549/*Node numerics:*/
550/*FUNCTION Node::ApplyConstraints{{{1*/
551void Node::ApplyConstraint(int dof,double value){
552
553 int index;
554
555 /*Dof should be added in the s set, describing which
556 * dofs are constrained to a certain value (dirichlet boundary condition*/
557 DofInSSet(dof-1);
558 this->indexing.svalues[dof-1]=value;
559}
560/*}}}*/
561/*FUNCTION Node::RelaxConstraint{{{1*/
562void Node::RelaxConstraint(int dof){
563
564 /*Dof should be added to the f-set, and taken out of the s-set:*/
565 DofInFSet(dof-1);
566 this->indexing.svalues[dof-1]=NAN;
567}
568/*}}}*/
569/*FUNCTION Node::CreateVecSets {{{1*/
570void Node::CreateVecSets(Vec pv_g,Vec pv_f,Vec pv_s){
571
572 double gvalue=1.0; //all nodes are in the g set;
573 double value;
574
575 int i;
576
577 for(i=0;i<this->indexing.gsize;i++){
578
579 /*g set: */
580 VecSetValues(pv_g,1,&indexing.gdoflist[i],&gvalue,INSERT_VALUES);
581
582 /*f set: */
583 value=(double)this->indexing.f_set[i];
584 VecSetValues(pv_f,1,&indexing.gdoflist[i],&value,INSERT_VALUES);
585
586 /*s set: */
587 value=(double)this->indexing.s_set[i];
588 VecSetValues(pv_s,1,&indexing.gdoflist[i],&value,INSERT_VALUES);
589
590 }
591
592
593}
594/*}}}*/
595/*FUNCTION Node::CreateNodalConstraints{{{1*/
596void Node::CreateNodalConstraints(Vec ys){
597
598 int i;
599 double* values=NULL;
600 int count;
601
602 /*Recover values for s set and plug them in constraints vector: */
603 if(this->indexing.ssize){
604 values=(double*)xmalloc(this->indexing.ssize*sizeof(double));
605 count=0;
606 for(i=0;i<this->indexing.gsize;i++){
607 if(this->indexing.s_set[i]){
608 values[count]=this->indexing.svalues[i];
609 _assert_(!isnan(values[count]));
610 count++;
611 }
612 }
613
614 /*Add values into constraint vector: */
615 VecSetValues(ys,this->indexing.ssize,this->indexing.sdoflist,values,INSERT_VALUES);
616 }
617
618 /*Free ressources:*/
619 xfree((void**)&values);
620
621
622}
623/*}}}*/
624/*FUNCTION Node::DofInSSet {{{1*/
625void Node::DofInSSet(int dof){
626
627 /*Put dof for this node into the s set (ie, this dof will be constrained
628 * to a fixed value during computations. */
629
630 this->indexing.f_set[dof]=0; //n splits into f (for which we solve) and s (single point constraints)
631 this->indexing.s_set[dof]=1;
632}
633/*}}}*/
634/*FUNCTION Node::DofInFSet {{{1*/
635void Node::DofInFSet(int dof){
636
637 /*Put dof for this node into the f set (ie, this dof will NOT be constrained
638 * to a fixed value during computations. */
639
640 this->indexing.f_set[dof]=1;
641 this->indexing.s_set[dof]=0;
642}
643/*}}}*/
644/*FUNCTION Node::FreezeDof{{{1*/
645void Node::FreezeDof(int dof){
646
647 DofInSSet(dof-1); //with 0 displacement for this dof.
648
649}
650/*}}}*/
651/*FUNCTION Node::GetApproximation {{{1*/
652int Node::GetApproximation(){
653
654 int approximation;
655
656 /*recover parameters: */
657 inputs->GetInputValue(&approximation,ApproximationEnum);
658
659 return approximation;
660}
661/*}}}*/
662/*FUNCTION Node::GetConnectivity {{{1*/
663int Node::GetConnectivity(){
664
665 Vertex* vertex=NULL;
666 vertex=(Vertex*)hvertex->delivers();
667 return vertex->connectivity;
668}
669/*}}}*/
670/*FUNCTION Node::GetNumberOfDofs{{{1*/
671int Node::GetNumberOfDofs(int approximation_enum,int setenum){
672
673 /*Get number of degrees of freedom in a node, for a certain set (g,f or s-set)
674 *and for a certain approximation type: */
675
676 int i;
677 int numdofs=0;
678
679 if(approximation_enum==NoneApproximationEnum){
680 if (setenum==GsetEnum) numdofs=this->indexing.gsize;
681 else if (setenum==FsetEnum) numdofs=this->indexing.fsize;
682 else if (setenum==SsetEnum) numdofs=this->indexing.ssize;
683 else _error_("%s%s%s"," set of enum type ",EnumToStringx(setenum)," not supported yet!");
684 }
685 else{
686 if(setenum==GsetEnum){
687 if(this->indexing.doftype){
688 numdofs=0;
689 for(i=0;i<this->indexing.gsize;i++){
690 if(this->indexing.doftype[i]==approximation_enum) numdofs++;
691 }
692 }
693 else numdofs=this->indexing.gsize;
694 }
695 else if (setenum==FsetEnum){
696 if(this->indexing.doftype){
697 numdofs=0;
698 for(i=0;i<this->indexing.gsize;i++){
699 if((this->indexing.doftype[i]==approximation_enum) && (this->indexing.f_set[i])) numdofs++;
700 }
701 }
702 else numdofs=this->indexing.fsize;
703 }
704 else if (setenum==SsetEnum){
705 if(this->indexing.doftype){
706 numdofs=0;
707 for(i=0;i<this->indexing.gsize;i++){
708 if((this->indexing.doftype[i]==approximation_enum) && (this->indexing.s_set[i])) numdofs++;
709 }
710 }
711 else numdofs=this->indexing.ssize;
712 }
713 else _error_("%s%s%s"," set of enum type ",EnumToStringx(setenum)," not supported yet!");
714 }
715 return numdofs;
716}
717/*}}}*/
718/*FUNCTION Node::GetSigma {{{1*/
719double Node::GetSigma(){
720 Vertex* vertex=NULL;
721
722 vertex=(Vertex*)hvertex->delivers();
723 return vertex->sigma;
724}
725/*}}}*/
726/*FUNCTION Node::GetX {{{1*/
727double Node::GetX(){
728 Vertex* vertex=NULL;
729
730 vertex=(Vertex*)hvertex->delivers();
731 return vertex->x;
732}
733/*}}}*/
734/*FUNCTION Node::GetY {{{1*/
735double Node::GetY(){
736 Vertex* vertex=NULL;
737
738 vertex=(Vertex*)hvertex->delivers();
739 return vertex->y;
740}
741/*}}}*/
742/*FUNCTION Node::GetZ {{{1*/
743double Node::GetZ(){
744 Vertex* vertex=NULL;
745
746 vertex=(Vertex*)hvertex->delivers();
747 return vertex->z;
748}
749/*}}}*/
750/*FUNCTION Node::IsClone {{{1*/
751int Node::IsClone(){
752
753 return indexing.clone;
754
755}
756/*}}}*/
757/*FUNCTION Node::IsOnBed {{{1*/
758int Node::IsOnBed(){
759
760 bool onbed;
761
762 /*recover parameters: */
763 inputs->GetInputValue(&onbed,MeshVertexonbedEnum);
764
765 return onbed;
766}
767/*}}}*/
768/*FUNCTION Node::IsGrounded {{{1*/
769int Node::IsGrounded(){
770
771 bool onsheet;
772
773 /*recover parameters: */
774 inputs->GetInputValue(&onsheet,MaskVertexongroundediceEnum);
775
776 return onsheet;
777}
778/*}}}*/
779/*FUNCTION Node::IsFloating {{{1*/
780int Node::IsFloating(){
781
782 bool onshelf;
783
784 /*recover parameters: */
785 inputs->GetInputValue(&onshelf,MaskVertexonfloatingiceEnum);
786
787 return onshelf;
788}
789/*}}}*/
790/*FUNCTION Node::IsOnSurface {{{1*/
791int Node::IsOnSurface(){
792
793 bool onsurface;
794
795 /*recover parameters: */
796 inputs->GetInputValue(&onsurface,MeshVertexonsurfaceEnum);
797
798 return onsurface;
799}
800/*}}}*/
801/*FUNCTION Node::InputUpdateFromVector(double* vector, int name, int type){{{1*/
802void Node::InputUpdateFromVector(double* vector, int name, int type){
803
804 /*Nothing updated yet*/
805}
806/*}}}*/
807/*FUNCTION Node::InputUpdateFromVector(int* vector, int name, int type){{{1*/
808void Node::InputUpdateFromVector(int* vector, int name, int type){
809
810 /*Nothing updated yet*/
811}
812/*}}}*/
813/*FUNCTION Node::InputUpdateFromVector(bool* vector, int name, int type){{{1*/
814void Node::InputUpdateFromVector(bool* vector, int name, int type){
815
816 /*Nothing updated yet*/
817}
818/*}}}*/
819/*FUNCTION Node::InputUpdateFromVectorDakota(double* vector, int name, int type){{{1*/
820void Node::InputUpdateFromVectorDakota(double* vector, int name, int type){
821
822 /*Nothing updated yet*/
823}
824/*}}}*/
825/*FUNCTION Node::InputUpdateFromMatrixDakota(double* matrix, int nrows, int ncols, int name, int type){{{1*/
826void Node::InputUpdateFromMatrixDakota(double* matrix, int nrows, int ncols, int name, int type){
827
828 /*Nothing updated yet*/
829}
830/*}}}*/
831/*FUNCTION Node::InputUpdateFromVectorDakota(int* vector, int name, int type){{{1*/
832void Node::InputUpdateFromVectorDakota(int* vector, int name, int type){
833
834 /*Nothing updated yet*/
835}
836/*}}}*/
837/*FUNCTION Node::InputUpdateFromVectorDakota(bool* vector, int name, int type){{{1*/
838void Node::InputUpdateFromVectorDakota(bool* vector, int name, int type){
839
840 /*Nothing updated yet*/
841}
842/*}}}*/
843/*FUNCTION Node::InputUpdateFromConstant(double constant, int name){{{1*/
844void Node::InputUpdateFromConstant(double constant, int name){
845
846 /*Nothing updated yet*/
847}
848/*}}}*/
849/*FUNCTION Node::InputUpdateFromConstant(int constant, int name){{{1*/
850void Node::InputUpdateFromConstant(int constant, int name){
851
852 /*Nothing updated yet*/
853}
854/*}}}*/
855/*FUNCTION Node::InputUpdateFromConstant(bool constant, int name){{{1*/
856void Node::InputUpdateFromConstant(bool constant, int name){
857
858 /*Nothing updated yet*/
859}
860/*}}}*/
861/*FUNCTION Node::UpdateSpcs {{{1*/
862void Node::UpdateSpcs(double* ys){
863
864 int count=0;
865 int i;
866
867 count=0;
868 for(i=0;i<this->indexing.gsize;i++){
869 if(this->indexing.s_set[i]){
870 this->indexing.svalues[i]=ys[this->indexing.sdoflist[count]];
871 count++;
872 }
873 }
874}
875/*}}}*/
876/*FUNCTION Node::VecMerge {{{1*/
877void Node::VecMerge(Vec ug, double* vector_serial,int setenum){
878
879 double* values=NULL;
880 int* indices=NULL;
881 int count=0;
882 int i;
883
884 if(setenum==FsetEnum){
885 if(this->indexing.fsize){
886 indices=(int*)xmalloc(this->indexing.fsize*sizeof(int));
887 values=(double*)xmalloc(this->indexing.fsize*sizeof(double));
888
889 for(i=0;i<this->indexing.gsize;i++){
890 if(this->indexing.f_set[i]){
891 _assert_(vector_serial);
892 values[count]=vector_serial[this->indexing.fdoflist[count]];
893 indices[count]=this->indexing.gdoflist[i];
894 count++;
895 }
896 }
897
898 /*Add values into ug: */
899 VecSetValues(ug,this->indexing.fsize,indices,(const double*)values,INSERT_VALUES);
900 }
901 }
902 else if(setenum==SsetEnum){
903 if(this->indexing.ssize){
904 indices=(int*)xmalloc(this->indexing.ssize*sizeof(int));
905 values=(double*)xmalloc(this->indexing.ssize*sizeof(double));
906
907 for(i=0;i<this->indexing.gsize;i++){
908 if(this->indexing.s_set[i]){
909 _assert_(vector_serial);
910 values[count]=vector_serial[this->indexing.sdoflist[count]];
911 indices[count]=this->indexing.gdoflist[i];
912 count++;
913 }
914 }
915
916 /*Add values into ug: */
917 VecSetValues(ug,this->indexing.ssize,indices,(const double*)values,INSERT_VALUES);
918 }
919 }
920 else _error_("VecMerge can only merge from the s or f-set onto the g-set!");
921
922 /*Free ressources:*/
923 xfree((void**)&values);
924 xfree((void**)&indices);
925}
926/*}}}*/
927/*FUNCTION Node::VecReduce {{{1*/
928void Node::VecReduce(Vec vector, double* ug_serial,int setenum){
929
930 double* values=NULL;
931 int count=0;
932 int i;
933
934 if(setenum==FsetEnum){
935 if(this->indexing.fsize){
936 values=(double*)xmalloc(this->indexing.fsize*sizeof(double));
937
938 for(i=0;i<this->indexing.gsize;i++){
939 if(this->indexing.f_set[i]){
940 _assert_(ug_serial);
941 values[count]=ug_serial[this->indexing.gdoflist[i]];
942 count++;
943 }
944 }
945
946 /*Add values into ug: */
947 VecSetValues(vector,this->indexing.fsize,this->indexing.fdoflist,(const double*)values,INSERT_VALUES);
948 }
949 }
950 else if(setenum==SsetEnum){
951 if(this->indexing.ssize){
952 values=(double*)xmalloc(this->indexing.ssize*sizeof(double));
953
954 for(i=0;i<this->indexing.gsize;i++){
955 if(this->indexing.s_set[i]){
956 _assert_(ug_serial);
957 values[count]=ug_serial[this->indexing.gdoflist[i]];
958 count++;
959 }
960 }
961
962 /*Add values into ug: */
963 VecSetValues(vector,this->indexing.ssize,this->indexing.sdoflist,(const double*)values,INSERT_VALUES);
964 }
965 }
966 else _error_("VecReduce can only merge from the s or f-set onto the g-set!");
967
968 /*Free ressources:*/
969 xfree((void**)&values);
970}
971/*}}}*/
972
973/* DofObject routines:*/
974/*FUNCTION Node::DistributeDofs{{{1*/
975void Node::DistributeDofs(int* pdofcount,int setenum){
976
977 int i;
978 extern int my_rank;
979 int dofcount;
980
981 dofcount=*pdofcount;
982
983 /*Initialize: */
984 if(setenum==FsetEnum) this->indexing.InitSet(setenum);
985 if(setenum==SsetEnum) this->indexing.InitSet(setenum);
986
987
988 /*For clone nodfs, don't distribute dofs, we will get them from another cpu in UpdateCloneDofs!*/
989 if(indexing.clone){
990 return;
991 }
992
993 /*This node should distribute dofs for setenum set (eg, f_set or s_set), go ahead: */
994
995 if(setenum==GsetEnum){
996 for(i=0;i<this->indexing.gsize;i++){
997 indexing.gdoflist[i]=dofcount+i;
998 }
999 dofcount+=this->indexing.gsize;
1000 }
1001 else if(setenum==FsetEnum){
1002 for(i=0;i<this->indexing.fsize;i++){
1003 indexing.fdoflist[i]=dofcount+i;
1004 }
1005 dofcount+=this->indexing.fsize;
1006 }
1007 else if(setenum==SsetEnum){
1008 for(i=0;i<this->indexing.ssize;i++){
1009 indexing.sdoflist[i]=dofcount+i;
1010 }
1011 dofcount+=this->indexing.ssize;
1012 }
1013 else _error_("%s%s%s"," set of enum type ",EnumToStringx(setenum)," not supported yet!");
1014
1015
1016 /*Assign output pointers: */
1017 *pdofcount=dofcount;
1018
1019}
1020/*}}}*/
1021/*FUNCTION Node::Off_setDofs{{{1*/
1022void Node::OffsetDofs(int dofcount,int setenum){
1023
1024 int i;
1025 extern int my_rank;
1026
1027 if(indexing.clone){
1028 /*This node is a clone, don't off_set the dofs!: */
1029 return;
1030 }
1031
1032 /*This node should off_set the dofs, go ahead: */
1033 if(setenum==GsetEnum){
1034 for(i=0;i<this->indexing.gsize;i++) indexing.gdoflist[i]+=dofcount;
1035 }
1036 else if(setenum==FsetEnum){
1037 for(i=0;i<this->indexing.fsize;i++) indexing.fdoflist[i]+=dofcount;
1038 }
1039 else if(setenum==SsetEnum){
1040 for(i=0;i<this->indexing.ssize;i++) indexing.sdoflist[i]+=dofcount;
1041 }
1042 else _error_("%s%s%s"," set of enum type ",EnumToStringx(setenum)," not supported yet!");
1043}
1044/*}}}*/
1045/*FUNCTION Node::ShowTrueDofs{{{1*/
1046void Node::ShowTrueDofs(int* truedofs, int ncols,int setenum){
1047
1048 int j;
1049 extern int my_rank;
1050
1051 /*Are we a clone? : */
1052 if(indexing.clone)return;
1053
1054 /*Ok, we are not a clone, just plug our dofs into truedofs: */
1055 if(setenum==GsetEnum)for(j=0;j<this->indexing.gsize;j++) *(truedofs+ncols*sid+j)=indexing.gdoflist[j];
1056 else if(setenum==FsetEnum)for(j=0;j<this->indexing.fsize;j++) *(truedofs+ncols*sid+j)=indexing.fdoflist[j];
1057 else if(setenum==SsetEnum)for(j=0;j<this->indexing.ssize;j++) *(truedofs+ncols*sid+j)=indexing.sdoflist[j];
1058 else _error_("%s%s%s"," set of enum type ",EnumToStringx(setenum)," not supported yet!");
1059
1060}
1061/*}}}*/
1062/*FUNCTION Node::UpdateCloneDofs{{{1*/
1063void Node::UpdateCloneDofs(int* alltruedofs,int ncols,int setenum){
1064
1065 int j;
1066 extern int my_rank;
1067
1068 /*If we are not a clone, don't update, we already have dofs!: */
1069 if(indexing.clone==0)return;
1070
1071
1072 /*Ok, we are a clone node, but we did not create the dofs for this node.
1073 * * Therefore, our doflist is garbage right now. Go pick it up in the alltruedofs: */
1074 if(setenum==GsetEnum)for(j=0;j<this->indexing.gsize;j++) indexing.gdoflist[j]=*(alltruedofs+ncols*sid+j);
1075 else if(setenum==FsetEnum)for(j=0;j<this->indexing.fsize;j++) indexing.fdoflist[j]=*(alltruedofs+ncols*sid+j);
1076 else if(setenum==SsetEnum)for(j=0;j<this->indexing.ssize;j++) indexing.sdoflist[j]=*(alltruedofs+ncols*sid+j);
1077 else _error_("%s%s%s"," set of enum type ",EnumToStringx(setenum)," not supported yet!");
1078
1079}
1080/*}}}*/
1081/*FUNCTION Node::SetClone {{{1*/
1082void Node::SetClone(int* minranks){
1083
1084 extern int my_rank;
1085
1086 if (minranks[sid]==my_rank){
1087 indexing.clone=0;
1088 }
1089 else{
1090 /*!there is a cpu with lower rank that has the same node,
1091 therefore, I am a clone*/
1092 indexing.clone=1;
1093 }
1094
1095}
1096/*}}}*/
Note: See TracBrowser for help on using the repository browser.