source: issm/trunk/src/c/Container/Elements.cpp@ 9878

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

Added support for P0 for results_on_vertices

File size: 7.9 KB
Line 
1/*
2 * \file Elements.c
3 * \brief: implementation of the Elements class, derived from DataSet class
4 */
5
6/*Headers: {{{1*/
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 <vector>
14#include <functional>
15#include <algorithm>
16#include <iostream>
17
18#include "./DataSet.h"
19#include "../shared/shared.h"
20#include "../include/include.h"
21#include "../EnumDefinitions/EnumDefinitions.h"
22
23using namespace std;
24/*}}}*/
25
26/*Object constructors and destructor*/
27/*FUNCTION Elements::Elements(){{{1*/
28Elements::Elements(){
29 return;
30}
31/*}}}*/
32/*FUNCTION Elements::Elements(int in_enum){{{1*/
33Elements::Elements(int in_enum): DataSet(in_enum){
34 //do nothing;
35 return;
36}
37/*}}}*/
38/*FUNCTION Elements::~Elements(){{{1*/
39Elements::~Elements(){
40 return;
41}
42/*}}}*/
43
44/*Object management*/
45/*FUNCTION Elements::Configure{{{1*/
46void Elements::Configure(Elements* elements,Loads* loads, Nodes* nodes, Vertices* vertices, Materials* materials,Parameters* parameters){
47
48 vector<Object*>::iterator object;
49 Element* element=NULL;
50
51 for ( object=objects.begin() ; object < objects.end(); object++ ){
52
53 element=(Element*)(*object);
54 element->Configure(elements,loads,nodes,materials,parameters);
55
56 }
57
58}
59/*}}}*/
60/*FUNCTION Elements::ProcessResultsUnits{{{1*/
61void Elements::ProcessResultsUnits(void){
62
63 //Process results to be output in the correct units
64 for(int i=0;i<this->Size();i++){
65 Element* element=(Element*)this->GetObjectByOffset(i);
66 element->ProcessResultsUnits();
67 }
68}
69/*}}}*/
70/*FUNCTION Elements::DeleteResults{{{1*/
71void Elements::DeleteResults(void){
72
73 for (int i=0;i<this->Size();i++){
74 Element* element=(Element*)this->GetObjectByOffset(i);
75 element->DeleteResults();
76 }
77}
78/*}}}*/
79/*FUNCTION Elements::ResultsToPatch{{{1*/
80Patch* Elements::ResultsToPatch(void){
81
82 /*output: */
83 Patch* patch=NULL;
84
85 /*intermediary: */
86 int i;
87 int count;
88 int numrows;
89 int numvertices;
90 int numnodes;
91 int max_numvertices;
92 int max_numnodes;
93 int element_numvertices;
94 int element_numrows;
95 int element_numnodes;
96
97 /*We are going to extract from the results within the elements, the desired results , and create a table
98 * of patch information, that will hold, for each element that computed the result that
99 * we desire, the enum_type of the result, the step and time, the id of the element, the interpolation type, the vertices ids, and the values
100 * at the nodes (could be different from the vertices). This will be used for visualization purposes.
101 * For example, we could build the following patch table, for velocities:
102 *
103 1. on a Beam element, Vx, at step 1, time .5, element id 1, interpolation type P0 (constant), vertices ids 1 and 2, one constant value 4.5
104 VxEnum 1 .5 1 P0 1 2 4.5 NaN NaN
105 2. on a Tria element, Vz, at step 2, time .8, element id 2, interpolation type P1 (linear), vertices ids 1 3 and 4, with values at 3 nodes 4.5, 3.2, 2.5
106 VzEnum 2 .8 2 P1 1 3 4 4.5 3.2 2.5
107 * ... etc ...
108 *
109 * So what do we need to build the table: the maximum number of vertices included in the table,
110 * and the maximum number of nodal values, as well as the number of rows. Once we have that,
111 * we ask the elements to fill their own row in the table, by looping on the elememnts.
112 *
113 * We will use the Patch object, which will store all of the information needed, and will be able
114 * to output itself to disk on its own. See object/Patch.h for format of this object.*/
115
116 /*First, determine maximum number of vertices, nodes, and number of results: */
117 numrows=0;
118 numvertices=0;
119 numnodes=0;
120
121 for(i=0;i<this->Size();i++){
122
123 Element* element=(Element*)this->GetObjectByOffset(i);
124 element->PatchSize(&element_numrows,&element_numvertices,&element_numnodes);
125
126 numrows+=element_numrows;
127 if(element_numvertices>numvertices)numvertices=element_numvertices;
128 if(element_numnodes>numnodes)numnodes=element_numnodes;
129 }
130
131 #ifdef _PARALLEL_
132 /*Synchronize across cluster, so as to not end up with different sizes for each patch on each cpu: */
133 MPI_Reduce (&numvertices,&max_numvertices,1,MPI_INT,MPI_MAX,0,MPI_COMM_WORLD );
134 MPI_Bcast(&max_numvertices,1,MPI_INT,0,MPI_COMM_WORLD);
135 numvertices=max_numvertices;
136
137 MPI_Reduce (&numnodes,&max_numnodes,1,MPI_INT,MPI_MAX,0,MPI_COMM_WORLD );
138 MPI_Bcast(&max_numnodes,1,MPI_INT,0,MPI_COMM_WORLD);
139 numnodes=max_numnodes;
140 #endif
141
142 /*Ok, initialize Patch object: */
143 patch=new Patch(numrows,numvertices,numnodes);
144
145 /*Now, go through elements, and fill the Patch object: */
146 count=0;
147 for(i=0;i<this->Size();i++){
148 Element* element=(Element*)this->GetObjectByOffset(i);
149 element->PatchFill(&count,patch);
150 }
151
152 return patch;
153}
154/*}}}*/
155/*FUNCTION Elements::SetCurrentConfiguration{{{1*/
156void Elements::SetCurrentConfiguration(Elements* elements,Loads* loads, Nodes* nodes, Vertices* vertices, Materials* materials,Parameters* parameters){
157
158 vector<Object*>::iterator object;
159 Element* element=NULL;
160
161 for ( object=objects.begin() ; object < objects.end(); object++ ){
162
163 element=(Element*)(*object);
164 element->SetCurrentConfiguration(elements,loads,nodes,materials,parameters);
165
166 }
167
168}
169/*}}}*/
170/*FUNCTION Elements::ToResults{{{1*/
171void Elements::ToResults(Results* results,Parameters* parameters,int step, double time){
172
173 extern int my_rank;
174
175 Patch *patch = NULL;
176 int *resultsenums = NULL;
177 double *vector_serial= NULL;
178 Vec vector = NULL;
179 bool io_gather;
180 bool results_on_vertices;
181 int numberofvertices;
182 int numberofresults;
183
184 /*Recover parameters: */
185 parameters->FindParam(&io_gather,SettingsIoGatherEnum);
186 parameters->FindParam(&results_on_vertices,SettingsResultsOnVerticesEnum);
187 parameters->FindParam(&numberofvertices,MeshNumberofverticesEnum);
188
189 if(results_on_vertices){
190 /*No patch here, we prepare vectors*/
191
192 /*OK, see what the first element of this partition has in stock (this is common to all partitions)*/
193 Element* element=(Element*)this->GetObjectByOffset(0);
194 element->ListResultsEnums(&resultsenums,&numberofresults);
195
196 /*Loop over all results and get nodal vector*/
197 for(int i=0;i<numberofresults;i++){
198 printf("%s\n",EnumToStringx(resultsenums[i]));
199
200 /*Get vector for result number i*/
201 vector=NewVec(numberofvertices);
202 for(int j=0;j<this->Size();j++){
203 Element* element=(Element*)this->GetObjectByOffset(j);
204 element->GetVectorFromResults(vector,resultsenums[i]);
205 }
206 VecAssemblyBegin(vector);
207 VecAssemblyEnd(vector);
208
209 /*Serialize and add to results*/
210 VecToMPISerial(&vector_serial,vector);
211 if(my_rank==0){
212 results->AddObject(new DoubleVecExternalResult(results->Size()+1,resultsenums[i],vector_serial,numberofvertices,step,time));
213 }
214
215 /*clean up*/
216 VecFree(&vector);
217 xfree((void**)&vector_serial);
218 }
219 }
220 else{
221 /*create patch object out of all results in this dataset: */
222 patch=this->ResultsToPatch();
223
224 /*Gather onto master cpu 0, if needed: */
225#ifdef _PARALLEL_
226 if(io_gather)patch->Gather();
227#endif
228
229 /*create result object and add to results dataset:*/
230 if (patch->maxvertices && patch->maxnodes){
231 results->AddObject(new IntExternalResult(results->Size()+1,PatchVerticesEnum,patch->maxvertices,step,time));
232 results->AddObject(new IntExternalResult(results->Size()+1,PatchNodesEnum, patch->maxnodes,step,time));
233 results->AddObject(new DoubleMatExternalResult(results->Size()+1,PatchEnum,patch->values,patch->numrows,patch->numcols,step,time));
234 }
235 }
236
237 /*Free ressources:*/
238 xfree((void**)&resultsenums);
239 delete patch;
240}
241/*}}}*/
242/*FUNCTION Elements::NumberOfElements{{{1*/
243int Elements::NumberOfElements(void){
244
245 int local_nelem=0;
246 int numberofelements;
247
248 #ifdef _PARALLEL_
249 local_nelem=this->Size();
250 MPI_Allreduce ( (void*)&local_nelem,(void*)&numberofelements,1,MPI_INT,MPI_SUM,MPI_COMM_WORLD);
251 #else
252 numberofelements=this->Size();
253 #endif
254
255 return numberofelements;
256}
257/*}}}*/
258/*FUNCTION Elements::InputCopy{{{1*/
259void Elements::InputDuplicate(int input_enum,int output_enum){
260
261 for(int i=0;i<this->Size();i++){
262 Element* element=(Element*)this->GetObjectByOffset(i);
263 element->InputDuplicate(input_enum,output_enum);
264 }
265}
266/*}}}*/
Note: See TracBrowser for help on using the repository browser.