1 // Copyright (C) 2007-2012 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
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 \defgroup noncoincidentdec NonCoincidentDEC
43 \section overview Overview
45 \c NonCoincidentDEC enables nonconservative remapping of fields
46 between two parallel codes.
47 The computation is possible for 3D meshes and 2D meshes.
48 It is not available for 3D surfaces. The computation enables fast parallel localization, and is based on a point in element search, followed
49 by a field evaluation at the point location. Thus, it is typically
50 faster than the \ref interpkerneldec which gives a \ref conservativeremapping.
51 It is particularly true for the initialisation phase (synchronize)
52 which is very computationnaly intensive in \ref interpkerneldec.
54 In the present version, only fields lying on elements are considered.
55 The value is estimated by locating the barycenter of the target
56 side cell in a source cell and sending the value of this source cell
57 as the value of the target cell.
59 \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."
61 \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."
63 A typical use of NonCoincidentDEC encompasses two distinct phases :
64 - A setup phase during which the intersection volumes are computed and the communication structures are setup. This corresponds to calling the NonCoincidentDEC::synchronize() method.
65 - 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.
67 The following code excerpt illutrates a typical use of the NonCoincidentDEC class.
71 NonCoincidentDEC dec(groupA, groupB);
72 dec.attachLocalField(field);
74 if (groupA.containsMyRank())
76 else if (groupB.containsMyRank())
81 Computing the field on the receiving side can be expressed in terms
82 of a matrix-vector product : \f$ \phi_t=W.\phi_s\f$, with \f$ \phi_t
83 \f$ the field on the target side and \f$ \phi_s \f$ the field on
85 In the P0-P0 case, this matrix is a plain rectangular matrix with one
86 non-zero element per row (with value 1). For instance, in the above figure, the matrix is :
89 \begin{tabular}{|cccc|}
98 fvm_nodal_t* medmemMeshToFVMMesh(const MEDMEM::MESH* mesh)
100 // create an FVM structure from the paramesh structure
101 fvm_nodal_t * fvm_nodal = fvm_nodal_create(mesh->getName().c_str(),mesh->getMeshDimension());
104 int nbtypes = mesh->getNumberOfTypes(MED_EN::MED_CELL);
105 const MED_EN::medGeometryElement* types = mesh->getTypes(MED_EN::MED_CELL);
106 for (int itype=0; itype<nbtypes; itype++)
108 fvm_element_t fvm_type;
109 switch (types[itype])
111 case MED_EN::MED_TRIA3 :
112 fvm_type=FVM_FACE_TRIA;
114 case MED_EN::MED_QUAD4 :
115 fvm_type=FVM_FACE_QUAD;
117 case MED_EN::MED_TETRA4 :
118 fvm_type=FVM_CELL_TETRA;
120 case MED_EN::MED_HEXA8 :
121 fvm_type=FVM_CELL_HEXA;
124 throw MEDEXCEPTION(" MED type conversion to fvm is not handled yet.");
129 fvm_lnum_t nbelems = mesh->getNumberOfElements(MED_EN::MED_CELL, types[itype]);
130 fvm_lnum_t* conn = new fvm_lnum_t[nbelems*(types[itype]%100)];
131 const int* mesh_conn =mesh->getConnectivity(MED_EN::MED_FULL_INTERLACE,MED_EN::MED_NODAL, MED_EN::MED_CELL, types[itype]);
132 for (int i=0; i<nbelems*(types[itype]%100); i++)
133 conn[i]=mesh_conn[i];
135 if (types[itype]==MED_EN::MED_TRIA3)
137 for (int i=0; i<nbelems;i++)
140 conn[3*i]=mesh_conn[3*i+1];
145 if (types[itype]==MED_EN::MED_TETRA4)
147 for (int i=0; i<nbelems;i++)
150 conn[4*i]=mesh_conn[4*i+1];
154 fvm_nodal_append_by_transfer(fvm_nodal, nbelems, fvm_type,0,0,0,conn,0);
156 int nbnodes= mesh->getNumberOfNodes();
157 int spacedim=mesh->getSpaceDimension();
158 fvm_coord_t* coords = new fvm_coord_t[nbnodes*spacedim];
159 const double* mesh_coords=mesh->getCoordinates(MED_EN::MED_FULL_INTERLACE);
160 for (int i=0; i<nbnodes*spacedim; i++)
161 coords[i]=mesh_coords[i];
162 fvm_nodal_transfer_vertices(fvm_nodal,coords);
167 fvm_nodal_t* medmemSupportToFVMMesh(const MEDMEM::SUPPORT* support)
170 // create an FVM structure from the paramesh structure
171 fvm_nodal_t * fvm_nodal = fvm_nodal_create(support->getName().c_str(),1);
173 const MEDMEM::MESH* mesh= support->getMesh();
176 MED_EN::medEntityMesh entity = support->getEntity();
178 int nbtypes = support->getNumberOfTypes();
179 const MED_EN::medGeometryElement* types = support->getTypes();
181 const int* type_offset = support->getNumberIndex();
183 //browsing through all types
184 for (int itype=0; itype<nbtypes; itype++)
186 fvm_element_t fvm_type;
187 switch (types[itype])
189 case MED_EN::MED_TRIA3 :
190 fvm_type=FVM_FACE_TRIA;
192 case MED_EN::MED_QUAD4 :
193 fvm_type=FVM_FACE_QUAD;
195 case MED_EN::MED_TETRA4 :
196 fvm_type=FVM_CELL_TETRA;
198 case MED_EN::MED_HEXA8 :
199 fvm_type=FVM_CELL_HEXA;
202 throw MEDEXCEPTION(" MED type conversion to fvm is not handled yet.");
206 fvm_lnum_t nbelems = support->getNumberOfElements(types[itype]);
208 //for a partial support, defining the element numbers that are taken into
209 //account in the support
210 fvm_lnum_t* elem_numbers=0;
211 if (!support->isOnAllElements())
213 elem_numbers = const_cast<fvm_lnum_t*> (support->getNumber(types[itype]));
215 //creating work arrays to store list of elems for partial suports
218 fvm_lnum_t* temp = new int[nbelems];
219 for (int i=0; i< nbelems; i++)
220 temp[i] = elem_numbers [i]-ioffset;
221 ioffset+=type_offset[itype];
225 //retrieving original mesh connectivity
226 fvm_lnum_t* conn = const_cast<fvm_lnum_t*> (mesh->getConnectivity(MED_EN::MED_FULL_INTERLACE,MED_EN::MED_NODAL,entity, types[itype]));
228 // adding the elements to the FVM structure
229 fvm_nodal_append_by_transfer(fvm_nodal, nbelems, fvm_type,0,0,0,conn,elem_numbers);
231 //cleaning work arrays (for partial supports)
232 if (!support->isOnAllElements() && itype>0)
233 delete[] elem_numbers;
239 NonCoincidentDEC::NonCoincidentDEC()
244 \addtogroup noncoincidentdec
248 /*! Constructor of a non coincident \ref 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 = ParaMEDMEM::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 = ParaMEDMEM::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();