1 // Copyright (C) 2007-2013 CEA/DEN, EDF R&D
3 // This library is free software; you can redistribute it and/or
4 // modify it under the terms of the GNU Lesser General Public
5 // License as published by the Free Software Foundation; either
6 // version 2.1 of the License.
8 // This library is distributed in the hope that it will be useful,
9 // but WITHOUT ANY WARRANTY; without even the implied warranty of
10 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
11 // Lesser General Public License for more details.
13 // You should have received a copy of the GNU Lesser General Public
14 // License along with this library; if not, write to the Free Software
15 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
17 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
19 #ifndef __MEDSPLITTER_MESHCOLLECTION_H__
20 #define __MEDSPLITTER_MESHCOLLECTION_H__
22 /*! Projects an integer or double field onto a new mesh
24 * \param old_collection initial mesh collection supporting the field
25 * \param fieldname name of the field to be cast
26 * \param itnumber time iteration number of the field to be cast
27 * \param ordernumber inner iteration number of the field to be cast
28 * \param type_of_template dummy argument passed for separate template instanciation
32 void MESHCollection::castFields (const MESHCollection& old_collection,
33 const string& fieldname, int itnumber, int ordernumber)
36 strcpy(fieldchar,fieldname.c_str());
38 int nb_old_domains=old_collection._topology->nbDomain();
39 int nb_new_domains=_topology->nbDomain();
41 vector <list<int> > element_array (_topology->nbDomain());
43 vector <MEDMEM::FIELD<T> *> old_fields ;
45 //cout << "MEDSPLITTER - reading fields from old collection"<<endl;
46 old_collection.getDriver()->readFields(old_fields,fieldchar, itnumber, ordernumber);
48 //if (dynamic_cast<MEDMEM::FIELD<int>*> (old_fields[0])==0)
49 // old_collection.getDriver()->readFieldsDouble(old_fields,fieldchar, itnumber, ordernumber);
51 // old_collection.getDriver()->readFieldsInt(old_fields,fieldchar, itnumber, ordernumber);
52 //cout <<"MEDSPLITTER - end of read"<<endl;
54 vector <const MEDMEM::SUPPORT*> old_supports(nb_old_domains);
55 vector <MEDMEM::SUPPORT*> new_supports(nb_new_domains);
56 vector <MEDMEM::FIELD<T> *> new_fields(nb_new_domains);
58 for (int iold = 0; iold < nb_old_domains; iold++)
60 old_supports[iold]=old_fields[iold]->getSupport();
62 for (int inew = 0; inew < nb_new_domains; inew++)
64 new_supports[inew]=new MEDMEM::SUPPORT();
67 //cout << "MEDSPLITTER - casting supports"<<endl;
68 castSupport(old_collection,old_supports,new_supports);
69 //cout << "MEDSPLITTER - end of cast"<<endl;
71 int nb_components = old_fields[0]->getNumberOfComponents();
72 const string* components_names = old_fields[0]->getComponentsNames();
73 const string* components_description = old_fields[0]->getComponentsDescriptions();
74 const string* components_units = old_fields[0]->getMEDComponentsUnits();
75 if (itnumber != old_fields[0]->getIterationNumber()) {cout << "PB with iteration number"<<endl;exit (1);}
76 int iteration_number=old_fields[0]->getIterationNumber();
77 int order_number=old_fields[0]->getOrderNumber();
78 double time=old_fields[0]->getTime();
79 bool has_gauss_pts = old_fields[0]->getGaussPresence();
80 //bool has_gauss_pts=true;
82 // Creating the fields objects
83 // Two different procedures are used whether the field contains Gauss points
86 vector <MEDMEM::MEDMEM_Array<T, FullInterlaceGaussPolicy>*> medarray (nb_new_domains);
88 for (int inew=0; inew < nb_new_domains; inew++)
91 new_fields[inew] = new MEDMEM::FIELD<T>(new_supports[inew],nb_components);
92 new_fields[inew]->setName(fieldname);
93 new_fields[inew]->setComponentsNames(components_names);
94 new_fields[inew]->setComponentsDescriptions(components_description);
95 new_fields[inew]->setMEDComponentsUnits(components_units);
96 new_fields[inew]->setIterationNumber(iteration_number);
97 new_fields[inew]->setOrderNumber(order_number);
98 new_fields[inew]->setTime(time);
102 new_fields[inew]=new MEDMEM::FIELD<T>();
103 //copying the structures describing the field
104 new_fields[inew]->setNumberOfComponents(nb_components);
105 new_fields[inew]->setSupport(new_supports[inew]);
106 new_fields[inew]->setName(fieldname);
107 new_fields[inew]->setComponentsNames(components_names);
108 new_fields[inew]->setComponentsDescriptions(components_description);
109 new_fields[inew]->setMEDComponentsUnits(components_units);
110 new_fields[inew]->setIterationNumber(iteration_number);
111 new_fields[inew]->setOrderNumber(order_number);
112 new_fields[inew]->setTime(time);
114 //counters for the gauss points
115 //nbtypegeo is the number of geometric types on the field
116 //nbelgeoc is the count of element for each type
117 //nbgaussgeo is the number of gauss points for each type
122 vector<int> gauss_pts_number;
125 // the GaussLocalization structures are browsed in the old field
126 // and copied to the new one
127 // the nbtypegeo counter is incremented so that
128 // it contains the number of types for which a gauss localization is defined
129 MED_EN::MESH_ENTITIES::const_iterator currentEntity;
130 std::list<MED_EN::medGeometryElement>::const_iterator iter;
131 currentEntity = MED_EN::meshEntities.find(MED_EN::MED_CELL);
132 for (iter = (*currentEntity).second.begin();iter != (*currentEntity).second.end(); iter++)
134 // getGaussLocalization throws an exception
135 // if the GaussLocalization is not available
138 MEDMEM::GAUSS_LOCALIZATION<FullInterlace>
139 gaussloc(old_fields[0]->getGaussLocalization(*iter));
140 new_fields[inew]->setGaussLocalization(*iter,gaussloc);
143 gauss_pts_number.push_back(gaussloc.getNbGauss());
145 nb_elem.push_back(new_supports[inew]->getNumberOfElements(*iter));
153 // les tableaux nbelgeoc commencent a 1
154 nbelgeoc = new int [nbtypegeo+1];
155 nbgaussgeo= new int [nbtypegeo+1];
159 for (int i=1; i<=nbtypegeo;i++)
162 nbelgeoc [i]=nb_elem[i-1]+nbelgeoc[i-1];
163 nbgaussgeo[i]=gauss_pts_number[i-1];
166 //a MEDMEM_Array structure is created to contain
167 //the data on the Gauss points
169 medarray[inew]=new MEDMEM::MEDMEM_Array<T,FullInterlaceGaussPolicy>
170 (new_fields[inew]->getNumberOfComponents(), size, nbtypegeo,
171 static_cast<const int* const>(nbelgeoc),
172 static_cast<const int* const> (nbgaussgeo));
173 new_fields[inew]->setArray(medarray[inew]);
174 // delete[] nbelgeoc;
175 // delete[] nbgaussgeo;
179 for (int idomain=0; idomain < old_collection._topology->nbDomain(); idomain++)
181 //retrieves the group igroup on domain idomain
183 const MEDMEM::SUPPORT* support = old_supports[idomain];
184 int nbelem = support->getNumberOfElements(MED_EN::MED_ALL_ELEMENTS);
185 if (nbelem==0) continue;
186 int* list_of_elems=0;
188 if (support->isOnAllElements())
190 list_of_elems = new int[nbelem];
191 for (int i=0; i<nbelem;i++)
192 list_of_elems[i]=i+1;
195 list_of_elems = const_cast<int*> (support->getNumber(MED_EN::MED_ALL_ELEMENTS));
197 int* array=new int[nbelem];
200 int* initial_local=0;
204 MED_EN::medEntityMesh entity = support->getEntity();
207 case MED_EN::MED_CELL :
209 local=new int[nbelem];
210 initial_local=list_of_elems;
211 old_collection.getTopology()->convertCellToGlobal(idomain,list_of_elems,nbelem,array);
212 _topology->convertGlobalCellList(array,nbelem,local,ip);
215 case MED_EN::MED_FACE :
216 old_collection.getTopology()->convertFaceToGlobal(idomain,list_of_elems,nbelem,array);
217 _topology->convertGlobalFaceListWithTwins(array,nbelem,local,ip,full_array,size);
218 initial_local=new int[size];
219 old_collection.getTopology()->convertGlobalFaceList(full_array,size,initial_local,idomain);
222 case MED_EN::MED_NODE :
223 old_collection.getTopology()->convertNodeToGlobal(idomain,list_of_elems,nbelem,array);
224 _topology->convertGlobalNodeListWithTwins(array,nbelem,local,ip,full_array,size);
225 initial_local=new int[size];
226 old_collection.getTopology()->convertGlobalNodeList(full_array,size,initial_local,idomain);
231 if (!has_gauss_pts) {
232 for (int i=0; i<size; i++)
234 for (int j=0; j<nb_components;j++)
236 T value = old_fields[idomain]->getValueIJ(initial_local[i],j+1);
238 new_fields[ip[i]]->setValueIJ(local[i],j+1, value);
244 for (int i=0; i<size; i++)
246 MED_EN::medGeometryElement type =
247 old_collection._mesh[idomain]->getElementType(entity,initial_local[i]);
248 int nb_gauss_points=old_fields[idomain]->getNumberOfGaussPoints(type);
249 for (int j=0; j<nb_components;j++)
250 for (int k=0;k<nb_gauss_points; k++)
252 T value = old_fields[idomain]->getValueIJK(initial_local[i],j+1,k+1);
253 medarray[ip[i]]->setIJK(local[i],j+1,k+1, value);
260 if (support->isOnAllElements())
262 delete[] list_of_elems;
265 if (entity==MED_EN::MED_FACE || entity==MED_EN::MED_NODE) delete[] initial_local;
268 retrieveDriver()->writeFields(new_fields,fieldchar);
269 //if (dynamic_cast<MEDMEM::FIELD<int>*>(new_fields[0])==0)
270 // retrieveDriver()->writeFieldsDouble(new_fields,fieldchar);
272 // retrieveDriver()->writeFieldsInt(new_fields,fieldchar);
273 for (int i=0; i<nb_new_domains; i++)
275 new_fields[i]->removeReference();
276 new_supports[i]->removeReference();
277 //delete medarray[i];
279 for (unsigned i=0; i<old_fields.size(); i++)
281 cout << "old field deletion" <<endl;
282 old_fields[i]->removeReference();
283 //old_supports[i]->removeReference();//as we did not addReference() neither created it
287 #endif /*MEDSPLITTER_MESHCOLLECTION_H_*/