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

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

merged trunk-jpl and trunk for revision 13393

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