1 // Copyright (C) 2007-2016 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, or (at your option) any later version.
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
21 #include "CommInterface.hxx"
22 #include "Topology.hxx"
23 #include "BlockTopology.hxx"
24 #include "ComponentTopology.hxx"
25 #include "ParaFIELD.hxx"
26 #include "MPIProcessorGroup.hxx"
28 #include "NonCoincidentDEC.hxx"
31 #include <fvm_parall.h>
32 #include <fvm_nodal.h>
33 #include <fvm_nodal_append.h>
34 #include <fvm_locator.h>
41 \anchor NonCoincidentDEC-det
42 \class NonCoincidentDEC
44 \c NonCoincidentDEC enables non-conservative remapping of fields
45 between two parallel codes.
46 The computation is possible for 3D meshes and 2D meshes.
47 It is not available for 3D surfaces.
49 The computation enables fast parallel localization, and is based on a point in element search, followed
50 by a field evaluation at the point location. Thus, it is typically
51 faster than the \ref InterpKernelDEC-det "InterpKernelDEC" which uses a
52 \ref InterpKerRemapGlobal "conservative remapping" (i.e. the same algorithms of volume
53 intersection as in the \ref remapper "sequential remapper")
54 It is particularly true for the initialisation phase (synchronize() method)
55 which has a significant computation cost in \ref InterpKernelDEC-det.
57 In the present version, only fields lying on elements ("P0") are considered.
58 The value is estimated by locating the barycenter of the target
59 side cell in a source cell and sending the value of this source cell
60 as the value of the target cell.
62 \image html NonCoincident_small.png "Example showing the transfer from a field based on a quadrangular mesh to a triangular mesh. The triangle barycenters are computed and located in the quadrangles. In a P0-P0 interpolation, the value on the quadrangle is then applied to the triangles whose barycenter lies within."
64 \image latex NonCoincident_small.eps "Example showing the transfer from a field based on a quadrangular mesh to a triangular mesh. The triangle barycenters are computed and located in the quadrangles. In a P0-P0 interpolation, the value on the quadrangle is then applied to the triangles whose barycenter lies within."
66 A typical use of NonCoincidentDEC encompasses two distinct phases :
67 - A setup phase during which the intersection volumes are computed and the communication structures are setup. This corresponds to calling the NonCoincidentDEC::synchronize() method.
68 - A use phase during which the remappings are actually performed. This corresponds to the calls to sendData() and recvData() which actually trigger the data exchange. The data exchange are synchronous in the current version of the library so that recvData() and sendData() calls must be synchronized on code A and code B processor groups.
70 The following code excerpt illutrates a typical use of the NonCoincidentDEC class.
74 NonCoincidentDEC dec(groupA, groupB);
75 dec.attachLocalField(field);
77 if (groupA.containsMyRank())
79 else if (groupB.containsMyRank())
84 Computing the field on the receiving side can be expressed in terms
85 of a matrix-vector product : \f$ \phi_t=W.\phi_s\f$, with \f$ \phi_t
86 \f$ the field on the target side and \f$ \phi_s \f$ the field on
88 In the P0-P0 case, this matrix is a plain rectangular matrix with one
89 non-zero element per row (with value 1). For instance, in the above figure, the matrix is :
92 \begin{tabular}{|cccc|}
101 fvm_nodal_t* medmemMeshToFVMMesh(const MEDMEM::MESH* mesh)
103 // create an FVM structure from the paramesh structure
104 std::string meshName(mesh->getName());//this line avoid that mesh->getName() object killed before fvm_nodal_create read the const char *.
105 fvm_nodal_t * fvm_nodal = fvm_nodal_create(meshName.c_str(),mesh->getMeshDimension());
108 int nbtypes = mesh->getNumberOfTypes(MED_EN::MED_CELL);
109 const MED_EN::medGeometryElement* types = mesh->getTypes(MED_EN::MED_CELL);
110 for (int itype=0; itype<nbtypes; itype++)
112 fvm_element_t fvm_type;
113 switch (types[itype])
115 case MED_EN::MED_TRIA3 :
116 fvm_type=FVM_FACE_TRIA;
118 case MED_EN::MED_QUAD4 :
119 fvm_type=FVM_FACE_QUAD;
121 case MED_EN::MED_TETRA4 :
122 fvm_type=FVM_CELL_TETRA;
124 case MED_EN::MED_HEXA8 :
125 fvm_type=FVM_CELL_HEXA;
128 throw MEDEXCEPTION(" MED type conversion to fvm is not handled yet.");
133 fvm_lnum_t nbelems = mesh->getNumberOfElements(MED_EN::MED_CELL, types[itype]);
134 fvm_lnum_t* conn = new fvm_lnum_t[nbelems*(types[itype]%100)];
135 const int* mesh_conn =mesh->getConnectivity(MED_EN::MED_FULL_INTERLACE,MED_EN::MED_NODAL, MED_EN::MED_CELL, types[itype]);
136 for (int i=0; i<nbelems*(types[itype]%100); i++)
137 conn[i]=mesh_conn[i];
139 if (types[itype]==MED_EN::MED_TRIA3)
141 for (int i=0; i<nbelems;i++)
144 conn[3*i]=mesh_conn[3*i+1];
149 if (types[itype]==MED_EN::MED_TETRA4)
151 for (int i=0; i<nbelems;i++)
154 conn[4*i]=mesh_conn[4*i+1];
158 fvm_nodal_append_by_transfer(fvm_nodal, nbelems, fvm_type,0,0,0,conn,0);
160 int nbnodes= mesh->getNumberOfNodes();
161 int spacedim=mesh->getSpaceDimension();
162 fvm_coord_t* coords = new fvm_coord_t[nbnodes*spacedim];
163 const double* mesh_coords=mesh->getCoordinates(MED_EN::MED_FULL_INTERLACE);
164 for (int i=0; i<nbnodes*spacedim; i++)
165 coords[i]=mesh_coords[i];
166 fvm_nodal_transfer_vertices(fvm_nodal,coords);
171 fvm_nodal_t* medmemSupportToFVMMesh(const MEDMEM::SUPPORT* support)
174 // create an FVM structure from the paramesh structure
175 std::string supportName(support->getName());//this line avoid that support->getName() object killed before fvm_nodal_create read the const char *.
176 fvm_nodal_t * fvm_nodal = fvm_nodal_create(supportName.c_str(),1);
178 const MEDMEM::MESH* mesh= support->getMesh();
181 MED_EN::medEntityMesh entity = support->getEntity();
183 int nbtypes = support->getNumberOfTypes();
184 const MED_EN::medGeometryElement* types = support->getTypes();
186 const int* type_offset = support->getNumberIndex();
188 //browsing through all types
189 for (int itype=0; itype<nbtypes; itype++)
191 fvm_element_t fvm_type;
192 switch (types[itype])
194 case MED_EN::MED_TRIA3 :
195 fvm_type=FVM_FACE_TRIA;
197 case MED_EN::MED_QUAD4 :
198 fvm_type=FVM_FACE_QUAD;
200 case MED_EN::MED_TETRA4 :
201 fvm_type=FVM_CELL_TETRA;
203 case MED_EN::MED_HEXA8 :
204 fvm_type=FVM_CELL_HEXA;
207 throw MEDEXCEPTION(" MED type conversion to fvm is not handled yet.");
211 fvm_lnum_t nbelems = support->getNumberOfElements(types[itype]);
213 //for a partial support, defining the element numbers that are taken into
214 //account in the support
215 fvm_lnum_t* elem_numbers=0;
216 if (!support->isOnAllElements())
218 elem_numbers = const_cast<fvm_lnum_t*> (support->getNumber(types[itype]));
220 //creating work arrays to store list of elems for partial suports
223 fvm_lnum_t* temp = new int[nbelems];
224 for (int i=0; i< nbelems; i++)
225 temp[i] = elem_numbers [i]-ioffset;
226 ioffset+=type_offset[itype];
230 //retrieving original mesh connectivity
231 fvm_lnum_t* conn = const_cast<fvm_lnum_t*> (mesh->getConnectivity(MED_EN::MED_FULL_INTERLACE,MED_EN::MED_NODAL,entity, types[itype]));
233 // adding the elements to the FVM structure
234 fvm_nodal_append_by_transfer(fvm_nodal, nbelems, fvm_type,0,0,0,conn,elem_numbers);
236 //cleaning work arrays (for partial supports)
237 if (!support->isOnAllElements() && itype>0)
238 delete[] elem_numbers;
244 NonCoincidentDEC::NonCoincidentDEC()
248 /*! Constructor of a non coincident \ref para-dec "DEC" with
249 * a source group on which lies a field lying on a mesh and a
250 * target group on which lies a mesh.
252 * \param source_group ProcessorGroup on the source side
253 * \param target_group ProcessorGroup on the target side
256 NonCoincidentDEC::NonCoincidentDEC(ProcessorGroup& source_group,
257 ProcessorGroup& target_group)
258 :DEC(source_group, target_group)
261 NonCoincidentDEC::~NonCoincidentDEC()
265 /*! Synchronization process. Calling this method
266 * synchronizes the topologies so that the target side
267 * gets the information which enable it to fetch the field value
268 * from the source side.
269 * A typical call is :
271 NonCoincidentDEC dec(source_group,target_group);
272 dec.attachLocalField(field);
276 void NonCoincidentDEC::synchronize()
279 //initializing FVM parallel environment
280 const MPI_Comm* comm=dynamic_cast<const MPIProcessorGroup*> (_union_group)->getComm();
281 fvm_parall_set_mpi_comm(*const_cast<MPI_Comm*> (comm));
284 //setting up the communication DEC on both sides
286 if (_source_group->containsMyRank())
288 MEDMEM::MESH* mesh = _local_field->getField()->getSupport()->getMesh();
289 fvm_nodal_t* source_nodal = MEDCoupling::medmemMeshToFVMMesh(mesh);
291 int target_size = _target_group->size() ;
292 int start_rank= _source_group->size();
293 const MPI_Comm* comm = (dynamic_cast<const MPIProcessorGroup*> (_union_group))->getComm();
295 _locator = fvm_locator_create(1e-6,
300 fvm_locator_set_nodal(_locator,
302 mesh->getSpaceDimension(),
308 _nb_distant_points = fvm_locator_get_n_dist_points(_locator);
309 _distant_coords = fvm_locator_get_dist_coords(_locator);
310 _distant_locations = fvm_locator_get_dist_locations(_locator);
313 if (_target_group->containsMyRank())
315 MEDMEM::MESH* mesh = _local_field->getField()->getSupport()->getMesh();
317 fvm_nodal_t* target_nodal = MEDCoupling::medmemMeshToFVMMesh(mesh);
318 int source_size = _source_group->size();
320 const MPI_Comm* comm = (dynamic_cast<const MPIProcessorGroup*> (_union_group))->getComm();
322 _locator = fvm_locator_create(1e-6,
326 int nbcells = mesh->getNumberOfElements(MED_EN::MED_CELL,MED_EN::MED_ALL_ELEMENTS);
327 const MEDMEM::SUPPORT* support=_local_field->getField()->getSupport();
328 MEDMEM::FIELD<double>* barycenter_coords = mesh->getBarycenter(support);
329 const double* coords = barycenter_coords->getValue();
330 fvm_locator_set_nodal(_locator,
332 mesh->getSpaceDimension(),
336 delete barycenter_coords;
341 /*! This method is called on the target group in order to
342 * trigger the retrieveal of field data. It must
343 * be called synchronously with a sendData() call on
346 void NonCoincidentDEC::recvData()
348 int nbelems = _local_field->getField()->getSupport()->getMesh()->getNumberOfElements(MED_EN::MED_CELL, MED_EN::MED_ALL_ELEMENTS);
349 int nbcomp = _local_field->getField()->getNumberOfComponents();
350 double* values = new double [nbelems*nbcomp];
351 fvm_locator_exchange_point_var(_locator,
358 _local_field->getField()->setValue(values);
359 if (_forced_renormalization_flag)
360 renormalizeTargetField();
364 /*! This method is called on the source group in order to
365 * send field data. It must be called synchronously with
366 * a recvData() call on
369 void NonCoincidentDEC::sendData()
371 const double* values=_local_field->getField()->getValue();
372 int nbcomp = _local_field->getField()->getNumberOfComponents();
373 double* distant_values = new double [_nb_distant_points*nbcomp];
375 //cheap interpolation : the value of the cell is transfered to the point
376 for (int i=0; i<_nb_distant_points; i++)
377 for (int j=0; j <nbcomp; j++)
378 distant_values[i*nbcomp+j]=values[(_distant_locations[i]-1)*nbcomp+j];
380 fvm_locator_exchange_point_var(_locator,
388 delete [] distant_values;
389 if (_forced_renormalization_flag)
390 renormalizeTargetField();