1 // Copyright (C) 2007-2015 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 nonconservative 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. The computation enables fast parallel localization, and is based on a point in element search, followed
48 by a field evaluation at the point location. Thus, it is typically
49 faster than the \ref InterpKernelDEC-det "InterpKernelDEC" which gives a
50 \ref InterpKerRemapGlobal "conservative remapping".
51 It is particularly true for the initialisation phase (synchronize)
52 which is very computationnaly intensive in \ref InterpKernelDEC-det.
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 std::string meshName(mesh->getName());//this line avoid that mesh->getName() object killed before fvm_nodal_create read the const char *.
102 fvm_nodal_t * fvm_nodal = fvm_nodal_create(meshName.c_str(),mesh->getMeshDimension());
105 int nbtypes = mesh->getNumberOfTypes(MED_EN::MED_CELL);
106 const MED_EN::medGeometryElement* types = mesh->getTypes(MED_EN::MED_CELL);
107 for (int itype=0; itype<nbtypes; itype++)
109 fvm_element_t fvm_type;
110 switch (types[itype])
112 case MED_EN::MED_TRIA3 :
113 fvm_type=FVM_FACE_TRIA;
115 case MED_EN::MED_QUAD4 :
116 fvm_type=FVM_FACE_QUAD;
118 case MED_EN::MED_TETRA4 :
119 fvm_type=FVM_CELL_TETRA;
121 case MED_EN::MED_HEXA8 :
122 fvm_type=FVM_CELL_HEXA;
125 throw MEDEXCEPTION(" MED type conversion to fvm is not handled yet.");
130 fvm_lnum_t nbelems = mesh->getNumberOfElements(MED_EN::MED_CELL, types[itype]);
131 fvm_lnum_t* conn = new fvm_lnum_t[nbelems*(types[itype]%100)];
132 const int* mesh_conn =mesh->getConnectivity(MED_EN::MED_FULL_INTERLACE,MED_EN::MED_NODAL, MED_EN::MED_CELL, types[itype]);
133 for (int i=0; i<nbelems*(types[itype]%100); i++)
134 conn[i]=mesh_conn[i];
136 if (types[itype]==MED_EN::MED_TRIA3)
138 for (int i=0; i<nbelems;i++)
141 conn[3*i]=mesh_conn[3*i+1];
146 if (types[itype]==MED_EN::MED_TETRA4)
148 for (int i=0; i<nbelems;i++)
151 conn[4*i]=mesh_conn[4*i+1];
155 fvm_nodal_append_by_transfer(fvm_nodal, nbelems, fvm_type,0,0,0,conn,0);
157 int nbnodes= mesh->getNumberOfNodes();
158 int spacedim=mesh->getSpaceDimension();
159 fvm_coord_t* coords = new fvm_coord_t[nbnodes*spacedim];
160 const double* mesh_coords=mesh->getCoordinates(MED_EN::MED_FULL_INTERLACE);
161 for (int i=0; i<nbnodes*spacedim; i++)
162 coords[i]=mesh_coords[i];
163 fvm_nodal_transfer_vertices(fvm_nodal,coords);
168 fvm_nodal_t* medmemSupportToFVMMesh(const MEDMEM::SUPPORT* support)
171 // create an FVM structure from the paramesh structure
172 std::string supportName(support->getName());//this line avoid that support->getName() object killed before fvm_nodal_create read the const char *.
173 fvm_nodal_t * fvm_nodal = fvm_nodal_create(supportName.c_str(),1);
175 const MEDMEM::MESH* mesh= support->getMesh();
178 MED_EN::medEntityMesh entity = support->getEntity();
180 int nbtypes = support->getNumberOfTypes();
181 const MED_EN::medGeometryElement* types = support->getTypes();
183 const int* type_offset = support->getNumberIndex();
185 //browsing through all types
186 for (int itype=0; itype<nbtypes; itype++)
188 fvm_element_t fvm_type;
189 switch (types[itype])
191 case MED_EN::MED_TRIA3 :
192 fvm_type=FVM_FACE_TRIA;
194 case MED_EN::MED_QUAD4 :
195 fvm_type=FVM_FACE_QUAD;
197 case MED_EN::MED_TETRA4 :
198 fvm_type=FVM_CELL_TETRA;
200 case MED_EN::MED_HEXA8 :
201 fvm_type=FVM_CELL_HEXA;
204 throw MEDEXCEPTION(" MED type conversion to fvm is not handled yet.");
208 fvm_lnum_t nbelems = support->getNumberOfElements(types[itype]);
210 //for a partial support, defining the element numbers that are taken into
211 //account in the support
212 fvm_lnum_t* elem_numbers=0;
213 if (!support->isOnAllElements())
215 elem_numbers = const_cast<fvm_lnum_t*> (support->getNumber(types[itype]));
217 //creating work arrays to store list of elems for partial suports
220 fvm_lnum_t* temp = new int[nbelems];
221 for (int i=0; i< nbelems; i++)
222 temp[i] = elem_numbers [i]-ioffset;
223 ioffset+=type_offset[itype];
227 //retrieving original mesh connectivity
228 fvm_lnum_t* conn = const_cast<fvm_lnum_t*> (mesh->getConnectivity(MED_EN::MED_FULL_INTERLACE,MED_EN::MED_NODAL,entity, types[itype]));
230 // adding the elements to the FVM structure
231 fvm_nodal_append_by_transfer(fvm_nodal, nbelems, fvm_type,0,0,0,conn,elem_numbers);
233 //cleaning work arrays (for partial supports)
234 if (!support->isOnAllElements() && itype>0)
235 delete[] elem_numbers;
241 NonCoincidentDEC::NonCoincidentDEC()
245 /*! Constructor of a non coincident \ref para-dec "DEC" with
246 * a source group on which lies a field lying on a mesh and a
247 * target group on which lies a mesh.
249 * \param source_group ProcessorGroup on the source side
250 * \param target_group ProcessorGroup on the target side
253 NonCoincidentDEC::NonCoincidentDEC(ProcessorGroup& source_group,
254 ProcessorGroup& target_group)
255 :DEC(source_group, target_group)
258 NonCoincidentDEC::~NonCoincidentDEC()
262 /*! Synchronization process. Calling this method
263 * synchronizes the topologies so that the target side
264 * gets the information which enable it to fetch the field value
265 * from the source side.
266 * A typical call is :
268 NonCoincidentDEC dec(source_group,target_group);
269 dec.attachLocalField(field);
273 void NonCoincidentDEC::synchronize()
276 //initializing FVM parallel environment
277 const MPI_Comm* comm=dynamic_cast<const MPIProcessorGroup*> (_union_group)->getComm();
278 fvm_parall_set_mpi_comm(*const_cast<MPI_Comm*> (comm));
281 //setting up the communication DEC on both sides
283 if (_source_group->containsMyRank())
285 MEDMEM::MESH* mesh = _local_field->getField()->getSupport()->getMesh();
286 fvm_nodal_t* source_nodal = ParaMEDMEM::medmemMeshToFVMMesh(mesh);
288 int target_size = _target_group->size() ;
289 int start_rank= _source_group->size();
290 const MPI_Comm* comm = (dynamic_cast<const MPIProcessorGroup*> (_union_group))->getComm();
292 _locator = fvm_locator_create(1e-6,
297 fvm_locator_set_nodal(_locator,
299 mesh->getSpaceDimension(),
305 _nb_distant_points = fvm_locator_get_n_dist_points(_locator);
306 _distant_coords = fvm_locator_get_dist_coords(_locator);
307 _distant_locations = fvm_locator_get_dist_locations(_locator);
310 if (_target_group->containsMyRank())
312 MEDMEM::MESH* mesh = _local_field->getField()->getSupport()->getMesh();
314 fvm_nodal_t* target_nodal = ParaMEDMEM::medmemMeshToFVMMesh(mesh);
315 int source_size = _source_group->size();
317 const MPI_Comm* comm = (dynamic_cast<const MPIProcessorGroup*> (_union_group))->getComm();
319 _locator = fvm_locator_create(1e-6,
323 int nbcells = mesh->getNumberOfElements(MED_EN::MED_CELL,MED_EN::MED_ALL_ELEMENTS);
324 const MEDMEM::SUPPORT* support=_local_field->getField()->getSupport();
325 MEDMEM::FIELD<double>* barycenter_coords = mesh->getBarycenter(support);
326 const double* coords = barycenter_coords->getValue();
327 fvm_locator_set_nodal(_locator,
329 mesh->getSpaceDimension(),
333 delete barycenter_coords;
338 /*! This method is called on the target group in order to
339 * trigger the retrieveal of field data. It must
340 * be called synchronously with a sendData() call on
343 void NonCoincidentDEC::recvData()
345 int nbelems = _local_field->getField()->getSupport()->getMesh()->getNumberOfElements(MED_EN::MED_CELL, MED_EN::MED_ALL_ELEMENTS);
346 int nbcomp = _local_field->getField()->getNumberOfComponents();
347 double* values = new double [nbelems*nbcomp];
348 fvm_locator_exchange_point_var(_locator,
355 _local_field->getField()->setValue(values);
356 if (_forced_renormalization_flag)
357 renormalizeTargetField();
361 /*! This method is called on the source group in order to
362 * send field data. It must be called synchronously with
363 * a recvData() call on
366 void NonCoincidentDEC::sendData()
368 const double* values=_local_field->getField()->getValue();
369 int nbcomp = _local_field->getField()->getNumberOfComponents();
370 double* distant_values = new double [_nb_distant_points*nbcomp];
372 //cheap interpolation : the value of the cell is transfered to the point
373 for (int i=0; i<_nb_distant_points; i++)
374 for (int j=0; j <nbcomp; j++)
375 distant_values[i*nbcomp+j]=values[(_distant_locations[i]-1)*nbcomp+j];
377 fvm_locator_exchange_point_var(_locator,
385 delete [] distant_values;
386 if (_forced_renormalization_flag)
387 renormalizeTargetField();