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

Last change on this file since 11684 was 11684, checked in by cborstad, 13 years ago

merged revisions 11428:11680 from trunk-jpl into branches/trunk-jpl-damage
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(Vector* 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 ys->SetValues(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(Vector* 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 ug->SetValues(this->indexing.fsize,indices,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 ug->SetValues(this->indexing.ssize,indices,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(Vector* 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 vector->SetValues(this->indexing.fsize,this->indexing.fdoflist,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 vector->SetValues(this->indexing.ssize,this->indexing.sdoflist,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.