source: issm/trunk/src/c/classes/Elements/Elements.cpp@ 15396

Last change on this file since 15396 was 15375, checked in by schlegel, 12 years ago

CHG: save z and thermal params so dakota will work with thermal
File size: 11.0 KB
Line 
1/*
2 * \file Elements.cpp
3 * \brief: Implementation of Elements class, derived from DataSet class
4 */
5
6/*Headers: {{{*/
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 "./Element.h"
14#include "./Elements.h"
15#include "../Params/Parameters.h"
16#include "../ExternalResults/Results.h"
17#include "../ExternalResults/GenericExternalResult.h"
18#include "../Patch.h"
19#include "../../toolkits/toolkits.h"
20#include "../../shared/shared.h"
21
22using namespace std;
23/*}}}*/
24
25/*Object constructors and destructor*/
26/*FUNCTION Elements::Elements(){{{*/
27Elements::Elements(){
28 enum_type=MeshElementsEnum;
29 return;
30}
31/*}}}*/
32/*FUNCTION Elements::~Elements(){{{*/
33Elements::~Elements(){
34 return;
35}
36/*}}}*/
37
38/*Object management*/
39/*FUNCTION Elements::Configure{{{*/
40void Elements::Configure(Elements* elements,Loads* loads, Nodes* nodes, Vertices* vertices, Materials* materials,Parameters* parameters){
41
42 vector<Object*>::iterator object;
43 Element* element=NULL;
44
45 for ( object=objects.begin() ; object < objects.end(); object++ ){
46
47 element=dynamic_cast<Element*>((*object));
48 element->Configure(elements,loads,nodes,vertices,materials,parameters);
49
50 }
51
52}
53/*}}}*/
54/*FUNCTION Elements::DeleteResults{{{*/
55void Elements::DeleteResults(void){
56
57 for (int i=0;i<this->Size();i++){
58 Element* element=dynamic_cast<Element*>(this->GetObjectByOffset(i));
59 element->DeleteResults();
60 }
61}
62/*}}}*/
63/*FUNCTION Elements::ResultsToPatch{{{*/
64Patch* Elements::ResultsToPatch(void){
65
66 /*output: */
67 Patch* patch=NULL;
68
69 /*intermediary: */
70 int i;
71 int count;
72 int numrows;
73 int numvertices;
74 int numnodes;
75 int max_numvertices;
76 int max_numnodes;
77 int element_numvertices;
78 int element_numrows;
79 int element_numnodes;
80
81 /*We are going to extract from the results within the elements, the desired results , and create a table
82 * of patch information, that will hold, for each element that computed the result that
83 * 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
84 * at the nodes (could be different from the vertices). This will be used for visualization purposes.
85 * For example, we could build the following patch table, for velocities:
86 *
87 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
88 VxEnum 1 .5 1 P0 1 2 4.5 NaN NaN
89 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
90 VzEnum 2 .8 2 P1 1 3 4 4.5 3.2 2.5
91 * ... etc ...
92 *
93 * So what do we need to build the table: the maximum number of vertices included in the table,
94 * and the maximum number of nodal values, as well as the number of rows. Once we have that,
95 * we ask the elements to fill their own row in the table, by looping on the elememnts.
96 *
97 * We will use the Patch object, which will store all of the information needed, and will be able
98 * to output itself to disk on its own. See object/Patch.h for format of this object.*/
99
100 /*First, determine maximum number of vertices, nodes, and number of results: */
101 numrows=0;
102 numvertices=0;
103 numnodes=0;
104
105 for(i=0;i<this->Size();i++){
106
107 Element* element=dynamic_cast<Element*>(this->GetObjectByOffset(i));
108 element->PatchSize(&element_numrows,&element_numvertices,&element_numnodes);
109
110 numrows+=element_numrows;
111 if(element_numvertices>numvertices)numvertices=element_numvertices;
112 if(element_numnodes>numnodes)numnodes=element_numnodes;
113 }
114
115 /*Synchronize across cluster, so as to not end up with different sizes for each patch on each cpu: */
116 #ifdef _HAVE_MPI_
117 MPI_Reduce (&numvertices,&max_numvertices,1,MPI_INT,MPI_MAX,0,IssmComm::GetComm() );
118 MPI_Bcast(&max_numvertices,1,MPI_INT,0,IssmComm::GetComm());
119 numvertices=max_numvertices;
120
121 MPI_Reduce (&numnodes,&max_numnodes,1,MPI_INT,MPI_MAX,0,IssmComm::GetComm() );
122 MPI_Bcast(&max_numnodes,1,MPI_INT,0,IssmComm::GetComm());
123 numnodes=max_numnodes;
124 #endif
125
126 /*Ok, initialize Patch object: */
127 patch=new Patch(numrows,numvertices,numnodes);
128
129 /*Now, go through elements, and fill the Patch object: */
130 count=0;
131 for(i=0;i<this->Size();i++){
132 Element* element=dynamic_cast<Element*>(this->GetObjectByOffset(i));
133 element->PatchFill(&count,patch);
134 }
135
136 return patch;
137}
138/*}}}*/
139/*FUNCTION Elements::SetCurrentConfiguration{{{*/
140void Elements::SetCurrentConfiguration(Elements* elements,Loads* loads, Nodes* nodes, Vertices* vertices, Materials* materials,Parameters* parameters){
141
142 vector<Object*>::iterator object;
143 Element* element=NULL;
144
145 for ( object=objects.begin() ; object < objects.end(); object++ ){
146
147 element=dynamic_cast<Element*>((*object));
148 element->SetCurrentConfiguration(elements,loads,nodes,materials,parameters);
149
150 }
151
152}
153/*}}}*/
154/*FUNCTION Elements::ToResults{{{*/
155void Elements::ToResults(Results* results,Parameters* parameters){
156
157 int my_rank;
158 int num_procs;
159
160 Patch *patch = NULL;
161 int *resultsenums = NULL;
162 int *resultssizes = NULL;
163 int *resultssteps = NULL;
164 IssmDouble *resultstimes = NULL;
165 IssmDouble *vector_serial = NULL;
166 Vector<IssmDouble> *vector = NULL;
167 bool io_gather;
168 bool results_as_patches;
169 int numberofvertices ,numberofelements;
170 int numberofresults ,vectorsize;
171 int rank;
172 int minrank;
173
174 /*recover my_rank:*/
175 my_rank=IssmComm::GetRank();
176 num_procs=IssmComm::GetSize();
177
178 /*Recover parameters: */
179 parameters->FindParam(&io_gather,SettingsIoGatherEnum);
180 parameters->FindParam(&results_as_patches,SettingsResultsAsPatchesEnum);
181 parameters->FindParam(&numberofvertices,MeshNumberofverticesEnum);
182 parameters->FindParam(&numberofelements,MeshNumberofelementsEnum);
183
184 if(!results_as_patches){
185 /*No patch here, we prepare vectors*/
186
187 /*Get rank of first cpu that has results*/
188 #ifdef _HAVE_MPI_
189 if(this->Size()) rank=my_rank;
190 else rank=num_procs;
191 MPI_Allreduce (&rank,&minrank,1,MPI_INT,MPI_MIN,IssmComm::GetComm());
192 #else
193 minrank=my_rank;
194 #endif
195
196 /*see what the first element of this partition has in stock (this is common to all partitions)*/
197 if(my_rank==minrank){
198 if(this->Size()==0) _error_("Cannot write results because there is no element??");
199 Element* element=dynamic_cast<Element*>(this->GetObjectByOffset(0));
200 element->ListResultsInfo(&resultsenums,&resultssizes,&resultstimes,&resultssteps,&numberofresults);
201 }
202 #ifdef _HAVE_MPI_
203 MPI_Bcast(&numberofresults,1,MPI_DOUBLE,minrank,IssmComm::GetComm());
204 #endif
205
206 /*Get out if there is no results. Otherwise broadcast info*/
207 if(!numberofresults) return;
208 #ifdef _HAVE_MPI_
209 if(my_rank!=minrank){
210 resultsenums=xNew<int>(numberofresults);
211 resultssizes=xNew<int>(numberofresults);
212 resultstimes=xNew<IssmDouble>(numberofresults);
213 resultssteps=xNew<int>(numberofresults);
214 }
215 MPI_Bcast(resultsenums,numberofresults,MPI_INT,minrank,IssmComm::GetComm());
216 MPI_Bcast(resultssizes,numberofresults,MPI_INT,minrank,IssmComm::GetComm());
217 MPI_Bcast(resultstimes,numberofresults,MPI_DOUBLE,minrank,IssmComm::GetComm());
218 MPI_Bcast(resultssteps,numberofresults,MPI_INT,minrank,IssmComm::GetComm());
219 #endif
220
221 /*Loop over all results and get nodal vector*/
222 for(int i=0;i<numberofresults;i++){
223
224 /*Get vector for result number i*/
225 if(resultssizes[i]==P1Enum) vectorsize=numberofvertices;
226 else if(resultssizes[i]==P0Enum) vectorsize=numberofelements;
227 else _error_("Unkown result size: " << EnumToStringx(resultssizes[i]));
228 vector=new Vector<IssmDouble>(vectorsize);
229
230 for(int j=0;j<this->Size();j++){
231 Element* element=dynamic_cast<Element*>(this->GetObjectByOffset(j));
232 element->GetVectorFromResults(vector,i,resultsenums[i],resultssizes[i]);
233 }
234 vector->Assemble();
235
236 /*Serialize and add to results*/
237 vector_serial=vector->ToMPISerial();
238 if(my_rank==0){
239 /*No need to add this vector for all cpus*/
240 #ifdef _HAVE_ADOLC_
241 IssmPDouble* vector_serial_passive=xNew<IssmPDouble>(vectorsize);
242 for(int k=0;k<vectorsize;k++)vector_serial_passive[k]=reCast<IssmPDouble>(vector_serial[k]);
243 results->AddObject(new GenericExternalResult<double*>(results->Size()+1,resultsenums[i],vector_serial_passive,vectorsize,1,resultssteps[i],resultstimes[i]));
244 xDelete<IssmPDouble>(vector_serial_passive);
245 #else
246 results->AddObject(new GenericExternalResult<double*>(results->Size()+1,resultsenums[i],vector_serial,vectorsize,1,resultssteps[i],resultstimes[i]));
247 #endif
248 }
249
250 /*clean up*/
251 delete vector;
252 xDelete<IssmDouble>(vector_serial);
253 }
254 }
255 else{
256 /*create patch object out of all results in this dataset: */
257 patch=this->ResultsToPatch();
258
259 /*Gather onto master cpu 0, if needed: */
260 if(io_gather)patch->Gather();
261
262 /*create result object and add to results dataset:*/
263 if (patch->maxvertices && patch->maxnodes){
264 results->AddObject(new GenericExternalResult<int>(results->Size()+1,PatchVerticesEnum,patch->maxvertices,1,0));
265 results->AddObject(new GenericExternalResult<int>(results->Size()+1,PatchNodesEnum, patch->maxnodes,1,0));
266 #ifdef _HAVE_ADOLC_
267 IssmPDouble* values_passive=xNew<IssmPDouble>(patch->numrows*patch->numcols);
268 for(int k=0;k<(patch->numrows*patch->numcols);k++)values_passive[k]=reCast<IssmPDouble>(patch->values[k]);
269 results->AddObject(new GenericExternalResult<double*>(results->Size()+1,PatchEnum, values_passive,patch->numrows,patch->numcols,1,0));
270 xDelete<IssmPDouble>(values_passive);
271 #else
272 results->AddObject(new GenericExternalResult<double*>(results->Size()+1,PatchEnum, patch->values,patch->numrows,patch->numcols,1,0));
273 #endif
274 }
275 }
276
277 /*Free ressources:*/
278 xDelete<int>(resultsenums);
279 xDelete<int>(resultssizes);
280 xDelete<int>(resultssteps);
281 xDelete<IssmDouble>(resultstimes);
282 delete patch;
283}
284/*}}}*/
285/*FUNCTION Elements::MaxNumNodes{{{*/
286int Elements::MaxNumNodes(void){
287
288 int max=0;
289 int allmax;
290 int numnodes=0;
291
292 /*Now go through all elements, and get how many nodes they own, unless they are clone nodes: */
293 for(int i=0;i<this->Size();i++){
294
295 Element* element=dynamic_cast<Element*>(this->GetObjectByOffset(i));
296 numnodes=element->GetNumberOfNodes();
297 if(numnodes>max)max=numnodes;
298 }
299
300 /*Grab max of all cpus: */
301 #ifdef _HAVE_MPI_
302 MPI_Allreduce((void*)&max,(void*)&allmax,1,MPI_INT,MPI_MAX,IssmComm::GetComm());
303 max=allmax;
304 #endif
305
306 return max;
307}
308/*}}}*/
309/*FUNCTION Elements::NumberOfElements{{{*/
310int Elements::NumberOfElements(void){
311
312 int local_nelem;
313 int numberofelements;
314
315 local_nelem=this->Size();
316 #ifdef _HAVE_MPI_
317 MPI_Allreduce ( (void*)&local_nelem,(void*)&numberofelements,1,MPI_INT,MPI_SUM,IssmComm::GetComm());
318 #else
319 numberofelements=local_nelem;
320 #endif
321
322 return numberofelements;
323}
324/*}}}*/
325/*FUNCTION Elements::InputDuplicate{{{*/
326void Elements::InputDuplicate(int input_enum,int output_enum){
327
328 for(int i=0;i<this->Size();i++){
329 Element* element=dynamic_cast<Element*>(this->GetObjectByOffset(i));
330 element->InputDuplicate(input_enum,output_enum);
331 }
332}
333/*}}}*/
Note: See TracBrowser for help on using the repository browser.