]> SALOME platform Git repositories - tools/medcoupling.git/blob - src/ParaMEDMEM/InterpKernelDEC.cxx
Salome HOME
Copyright update 2020
[tools/medcoupling.git] / src / ParaMEDMEM / InterpKernelDEC.cxx
1 // Copyright (C) 2007-2020  CEA/DEN, EDF R&D
2 //
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.
7 //
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.
12 //
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
16 //
17 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
18 //
19
20 #include <mpi.h>
21 #include "CommInterface.hxx"
22 #include "Topology.hxx"
23 #include "BlockTopology.hxx"
24 #include "ComponentTopology.hxx"
25 #include "ParaFIELD.hxx"
26 #include "MPIProcessorGroup.hxx"
27 #include "ParaMESH.hxx"
28 #include "DEC.hxx"
29 #include "InterpolationMatrix.hxx"
30 #include "InterpKernelDEC.hxx"
31 #include "ElementLocator.hxx"
32
33 namespace MEDCoupling
34 {  
35
36   /*!
37     \anchor InterpKernelDEC-det
38     \class InterpKernelDEC
39
40     \section InterpKernelDEC-over Overview
41
42     The InterpKernelDEC enables the \ref InterpKerRemapGlobal "remapping" (or interpolation) of fields between
43     two parallel codes.
44
45     The projection
46     methodology is based on the algorithms of %INTERP_KERNEL, that is to say, they work in a similar fashion than
47     what the \ref remapper "sequential remapper" does. The following \ref discretization "projection methods"
48     are supported: P0->P0 (the most common case), P1->P0, P0->P1.
49
50     The computation is possible for 3D meshes, 2D meshes, and 3D-surface
51     meshes. Dimensions must be identical for code A and code B (for instance, though it could be
52     desirable, it is not yet possible to couple 3D surfaces with 2D surfaces).
53
54     The name "InterpKernelDEC" comes from the fact that this class uses exactly the same algorithms
55     as the sequential remapper. Both this class and the sequential
56     \ref MEDCoupling::MEDCouplingRemapper "MEDCouplingRemapper" are built on top of the %INTERP_KERNEL
57     algorithms (notably the computation of the intersection volumes).
58
59     Among the important properties inherited from the parent abstract class \ref DisjointDEC-det "DisjointDEC",
60     the two \ref MPIProcessorGroup-det "processor groups" (source and target) must have a void intersection.
61
62     \image html NonCoincident_small.png "Transfer of a field supported by a quadrangular mesh to a triangular mesh".
63
64     \image latex NonCoincident_small.eps "Transfer of a field supported by a quadrangular mesh to a triangular mesh"
65
66     In the figure above we see the transfer of a field based on a quadrangular mesh to a new field supported by
67     a triangular mesh. In a P0-P0 interpolation, to obtain the value on a triangle, the values on the
68     quadrangles are weighted by their intersection area and summed.
69
70     A typical use of InterpKernelDEC encompasses two distinct phases :
71     - A setup phase during which the intersection volumes are computed and the communication structures are
72     setup. This corresponds to calling the InterpKernelDEC::synchronize() method.
73     - A running phase during which the projections are actually performed. This corresponds to the calls to
74     sendData() and recvData() which actually trigger the data exchange. The data exchange are synchronous
75     in the current version of the library so that recvData() and sendData() calls must be synchronized
76     on code A and code B processor groups.
77
78     The following code excerpt illustrates a typical use of the InterpKernelDEC class.
79
80     \code
81     ...
82     InterpKernelDEC dec(groupA, groupB);
83     dec.attachLocalField(field);
84     dec.synchronize();
85     if (groupA.containsMyRank())
86     dec.recvData();
87     else if (groupB.containsMyRank())
88     dec.sendData();
89     ...
90     \endcode
91     A \ref InterpKerRemapGlobal "remapping" of the field from the source mesh to the target mesh is performed by
92     the function synchronise(), which computes the interpolation matrix.
93
94     Computing the field on the receiving side can be expressed in terms of a matrix-vector product :
95     \f$ \phi_t=W.\phi_s\f$, with \f$ \phi_t \f$ the field on the target side and \f$ \phi_s \f$ the field
96     on the source side.
97     When remapping a 3D surface to another 3D surface, a projection phase is necessary to match elements
98     from both sides. Care must be taken when defining this projection to obtain a
99     \ref InterpKerRemapGlobal "conservative remapping".
100
101     In the P0-P0 case, this matrix is a plain rectangular matrix with coefficients equal to the
102     intersection areas between triangle and quadrangles. For instance, in the above figure, the matrix
103     is :
104
105     \f[
106     \begin{tabular}{|cccc|}
107     0.72 & 0 & 0.2 & 0 \\
108     0.46 & 0 & 0.51 & 0.03\\
109     0.42 & 0.53 & 0 & 0.05\\
110     0 & 0 & 0.92 & 0.05 \\
111     \end{tabular}
112     \f]
113
114     \section InterpKernelDEC-options Options
115     On top of the usual \ref MEDCoupling::DECOptions "DEC options", the options supported by %InterpKernelDEC objects are
116     related to the underlying \ref InterpKerIntersectors "intersector class".
117     All the options available in the intersector objects are
118     available for the %InterpKernelDEC object. The various options available for  intersectors can
119     be reviewed in \ref InterpKerIntersectors.
120  
121     For instance :
122     \verbatim
123     InterpKernelDEC dec(source_group, target_group);
124     dec.attachLocalField(field);
125     dec.setDoRotate(false);
126     dec.setPrecision(1e-12);
127     dec.synchronize();
128     \endverbatim
129
130     \warning{  Options must be set before calling the synchronize method. }
131   */
132   
133   InterpKernelDEC::InterpKernelDEC():
134     DisjointDEC(),
135     _nb_distant_points(0), _distant_coords(0),
136     _distant_locations(0), _interpolation_matrix(0)
137   {  
138   }
139
140   /*!
141     This constructor creates an InterpKernelDEC which has \a source_group as a working side 
142     and  \a target_group as an idle side. All the processors will actually participate, but intersection computations will be performed on the working side during the \a synchronize() phase.
143     The constructor must be called synchronously on all processors of both processor groups.
144     The source group and target group MUST form a partition of all the procs within the communicator passed as 'world_comm'
145     when building the group.
146
147     \param source_group working side ProcessorGroup
148     \param target_group lazy side ProcessorGroup
149
150   */
151   InterpKernelDEC::InterpKernelDEC(ProcessorGroup& source_group, ProcessorGroup& target_group):
152     DisjointDEC(source_group, target_group),
153     _nb_distant_points(0), _distant_coords(0),
154     _distant_locations(0), _interpolation_matrix(0)
155   {
156
157   }
158
159
160   /*!
161    * Creates an InterpKernelDEC from a set of source procs IDs and target group IDs.
162    * The difference with the ctor using groups is that the set of procs might not cover entirely MPI_COMM_WORLD
163    * (a sub-communicator holding the union of source and target procs is recreated internally).
164    */
165   InterpKernelDEC::InterpKernelDEC(const std::set<int>& src_ids, const std::set<int>& trg_ids,
166                                    const MPI_Comm& world_comm):
167     DisjointDEC(src_ids,trg_ids,world_comm),
168     _nb_distant_points(0), _distant_coords(0),
169     _distant_locations(0), _interpolation_matrix(0)
170   {
171   }
172
173   InterpKernelDEC::~InterpKernelDEC()
174   {
175     if (_interpolation_matrix !=0)
176       delete _interpolation_matrix;
177   } 
178
179   /*! 
180     \brief Synchronization process for exchanging topologies.
181
182     This method prepares all the structures necessary for sending data from a processor group to the other. It uses the mesh
183     underlying the fields that have been set with attachLocalField method.
184     It works in four steps :
185     -# Bounding boxes are computed for each sub-domain,
186     -# The lazy side mesh parts that are likely to intersect the working side local processor are sent to the working side,
187     -# The working side calls the interpolation kernel to compute the intersection between local and imported mesh.
188     -# The lazy side is updated so that it knows the structure of the data that will be sent by
189     the working side during a \a sendData() call.
190
191   */
192   void InterpKernelDEC::synchronize()
193   {
194     if(!isInUnion())
195       return ;
196     delete _interpolation_matrix;
197     _interpolation_matrix = new InterpolationMatrix (_local_field, *_source_group,*_target_group,*this,*this); 
198
199     //setting up the communication DEC on both sides  
200     if (_source_group->containsMyRank())
201       {
202         //locate the distant meshes
203         ElementLocator locator(*_local_field, *_target_group, *_source_group);
204         //transferring option from InterpKernelDEC to ElementLocator   
205         locator.copyOptions(*this);
206         MEDCouplingPointSet* distant_mesh=0; 
207         mcIdType* distant_ids=0;
208         std::string distantMeth;
209         for (int i=0; i<_target_group->size(); i++)
210           {
211             //        int idistant_proc = (i+_source_group->myRank())%_target_group->size();
212             int idistant_proc=i;
213
214             //gathers pieces of the target meshes that can intersect the local mesh
215             locator.exchangeMesh(idistant_proc,distant_mesh,distant_ids);
216             if (distant_mesh !=0)
217               {
218                 locator.exchangeMethod(_method,idistant_proc,distantMeth);
219                 //adds the contribution of the distant mesh on the local one
220                 int idistant_proc_in_union=_union_group->translateRank(_target_group,idistant_proc);
221                 //std::cout <<"add contribution from proc "<<idistant_proc_in_union<<" to proc "<<_union_group->myRank()<<std::endl;
222                 _interpolation_matrix->addContribution(*distant_mesh,idistant_proc_in_union,distant_ids,_method,distantMeth);
223                 distant_mesh->decrRef();
224                 delete [] distant_ids;
225                 distant_mesh=0;
226                 distant_ids=0;
227               }
228           }
229        _interpolation_matrix->finishContributionW(locator);
230       }
231
232     if (_target_group->containsMyRank())
233       {
234         ElementLocator locator(*_local_field, *_source_group, *_target_group);
235         //transferring option from InterpKernelDEC to ElementLocator
236         locator.copyOptions(*this);
237         MEDCouplingPointSet* distant_mesh=0;
238         mcIdType* distant_ids=0;
239         for (int i=0; i<_source_group->size(); i++)
240           {
241             //        int idistant_proc = (i+_target_group->myRank())%_source_group->size();
242             int  idistant_proc=i;
243             //gathers pieces of the target meshes that can intersect the local mesh
244             locator.exchangeMesh(idistant_proc,distant_mesh,distant_ids);
245             //std::cout << " Data sent from "<<_union_group->myRank()<<" to source proc "<< idistant_proc<<std::endl;
246             if (distant_mesh!=0)
247               {
248                 std::string distantMeth;
249                 locator.exchangeMethod(_method,idistant_proc,distantMeth);
250                 distant_mesh->decrRef();
251                 delete [] distant_ids;
252                 distant_mesh=0;
253                 distant_ids=0;
254               }
255           }
256         _interpolation_matrix->finishContributionL(locator);
257       }
258     _interpolation_matrix->prepare();
259   }
260
261
262   /*!
263     Receives the data whether the processor is on the working side or on the lazy side. It must match a \a sendData() call on the other side.
264   */
265   void InterpKernelDEC::recvData()
266   {
267     if (_source_group->containsMyRank())
268       _interpolation_matrix->transposeMultiply(*_local_field->getField());
269     else if (_target_group->containsMyRank())
270       {
271         _interpolation_matrix->multiply(*_local_field->getField());
272         if (getForcedRenormalization())
273           renormalizeTargetField(getMeasureAbsStatus());
274       }
275   }
276
277
278   /*!
279     Receives the data at time \a time in asynchronous mode. The value of the field
280     will be time-interpolated from the field values received.
281     \param time time at which the value is desired
282   */
283   void InterpKernelDEC::recvData( double time )
284   {
285     _interpolation_matrix->getAccessDEC()->setTime(time);
286     recvData() ;
287   }
288
289   /*!
290     Sends the data whether the processor is on the working side or on the lazy side.
291     It must match a recvData() call on the other side.
292   */
293   void InterpKernelDEC::sendData()
294   {
295     if (_source_group->containsMyRank())
296       {
297     
298         _interpolation_matrix->multiply(*_local_field->getField());
299         if (getForcedRenormalization())
300           renormalizeTargetField(getMeasureAbsStatus());
301     
302       }
303     else if (_target_group->containsMyRank())
304       _interpolation_matrix->transposeMultiply(*_local_field->getField());
305   }
306
307   /*!
308     Sends the data available at time \a time in asynchronous mode. 
309     \param time time at which the value is available
310     \param deltatime time interval between the value presently sent and the next one. 
311   */
312   void InterpKernelDEC::sendData( double time , double deltatime )
313   {
314     _interpolation_matrix->getAccessDEC()->setTime(time,deltatime);
315     sendData() ;
316   }
317   
318 }