[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 |
|
---|
| 23 | using namespace std;
|
---|
| 24 | /*}}}*/
|
---|
| 25 |
|
---|
| 26 | /*Object constructors and destructor*/
|
---|
[12706] | 27 | /*FUNCTION Elements::Elements(){{{*/
|
---|
[4216] | 28 | Elements::Elements(){
|
---|
[10522] | 29 | enum_type=MeshElementsEnum;
|
---|
[4216] | 30 | return;
|
---|
| 31 | }
|
---|
| 32 | /*}}}*/
|
---|
[12706] | 33 | /*FUNCTION Elements::~Elements(){{{*/
|
---|
[4216] | 34 | Elements::~Elements(){
|
---|
| 35 | return;
|
---|
| 36 | }
|
---|
| 37 | /*}}}*/
|
---|
| 38 |
|
---|
| 39 | /*Object management*/
|
---|
[12706] | 40 | /*FUNCTION Elements::Configure{{{*/
|
---|
[6411] | 41 | void 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] | 56 | void 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] | 66 | void 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] | 75 | Patch* 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] | 151 | void 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] | 166 | void 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] | 295 | int 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] | 311 | void 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 | /*}}}*/
|
---|