]> SALOME platform Git repositories - tools/medcoupling.git/blob - src/ParaMEDMEM/InterpKernelDEC.cxx
Salome HOME
Get relevant changes from V7_dev branch (copyright update, adm files etc)
[tools/medcoupling.git] / src / ParaMEDMEM / InterpKernelDEC.cxx
1 // Copyright (C) 2007-2016  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
145     \param source_group working side ProcessorGroup
146     \param target_group lazy side ProcessorGroup
147
148   */
149   InterpKernelDEC::InterpKernelDEC(ProcessorGroup& source_group, ProcessorGroup& target_group):
150     DisjointDEC(source_group, target_group),
151     _nb_distant_points(0), _distant_coords(0),
152     _distant_locations(0), _interpolation_matrix(0)
153   {
154
155   }
156
157   InterpKernelDEC::InterpKernelDEC(const std::set<int>& src_ids, const std::set<int>& trg_ids,
158                                    const MPI_Comm& world_comm):
159     DisjointDEC(src_ids,trg_ids,world_comm),
160     _nb_distant_points(0), _distant_coords(0),
161     _distant_locations(0), _interpolation_matrix(0)
162   {
163   }
164
165   InterpKernelDEC::~InterpKernelDEC()
166   {
167     if (_interpolation_matrix !=0)
168       delete _interpolation_matrix;
169   } 
170
171   /*! 
172     \brief Synchronization process for exchanging topologies.
173
174     This method prepares all the structures necessary for sending data from a processor group to the other. It uses the mesh
175     underlying the fields that have been set with attachLocalField method.
176     It works in four steps :
177     -# Bounding boxes are computed for each sub-domain,
178     -# The lazy side mesh parts that are likely to intersect the working side local processor are sent to the working side,
179     -# The working side calls the interpolation kernel to compute the intersection between local and imported mesh.
180     -# The lazy side is updated so that it knows the structure of the data that will be sent by
181     the working side during a \a sendData() call.
182
183   */
184   void InterpKernelDEC::synchronize()
185   {
186     if(!isInUnion())
187       return ;
188     delete _interpolation_matrix;
189     _interpolation_matrix = new InterpolationMatrix (_local_field, *_source_group,*_target_group,*this,*this); 
190
191     //setting up the communication DEC on both sides  
192     if (_source_group->containsMyRank())
193       {
194         //locate the distant meshes
195         ElementLocator locator(*_local_field, *_target_group, *_source_group);
196         //transfering option from InterpKernelDEC to ElementLocator   
197         locator.copyOptions(*this);
198         MEDCouplingPointSet* distant_mesh=0; 
199         int* distant_ids=0;
200         std::string distantMeth;
201         for (int i=0; i<_target_group->size(); i++)
202           {
203             //        int idistant_proc = (i+_source_group->myRank())%_target_group->size();
204             int idistant_proc=i;
205
206             //gathers pieces of the target meshes that can intersect the local mesh
207             locator.exchangeMesh(idistant_proc,distant_mesh,distant_ids);
208             if (distant_mesh !=0)
209               {
210                 locator.exchangeMethod(_method,idistant_proc,distantMeth);
211                 //adds the contribution of the distant mesh on the local one
212                 int idistant_proc_in_union=_union_group->translateRank(_target_group,idistant_proc);
213                 //std::cout <<"add contribution from proc "<<idistant_proc_in_union<<" to proc "<<_union_group->myRank()<<std::endl;
214                 _interpolation_matrix->addContribution(*distant_mesh,idistant_proc_in_union,distant_ids,_method,distantMeth);
215                 distant_mesh->decrRef();
216                 delete [] distant_ids;
217                 distant_mesh=0;
218                 distant_ids=0;
219               }
220           }
221        _interpolation_matrix->finishContributionW(locator);
222       }
223
224     if (_target_group->containsMyRank())
225       {
226         ElementLocator locator(*_local_field, *_source_group, *_target_group);
227         //transfering option from InterpKernelDEC to ElementLocator
228         locator.copyOptions(*this);
229         MEDCouplingPointSet* distant_mesh=0;
230         int* distant_ids=0;
231         for (int i=0; i<_source_group->size(); i++)
232           {
233             //        int idistant_proc = (i+_target_group->myRank())%_source_group->size();
234             int  idistant_proc=i;
235             //gathers pieces of the target meshes that can intersect the local mesh
236             locator.exchangeMesh(idistant_proc,distant_mesh,distant_ids);
237             //std::cout << " Data sent from "<<_union_group->myRank()<<" to source proc "<< idistant_proc<<std::endl;
238             if (distant_mesh!=0)
239               {
240                 std::string distantMeth;
241                 locator.exchangeMethod(_method,idistant_proc,distantMeth);
242                 distant_mesh->decrRef();
243                 delete [] distant_ids;
244                 distant_mesh=0;
245                 distant_ids=0;
246               }
247           }
248         _interpolation_matrix->finishContributionL(locator);
249       }
250     _interpolation_matrix->prepare();
251   }
252
253
254   /*!
255     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.
256   */
257   void InterpKernelDEC::recvData()
258   {
259     if (_source_group->containsMyRank())
260       _interpolation_matrix->transposeMultiply(*_local_field->getField());
261     else if (_target_group->containsMyRank())
262       {
263         _interpolation_matrix->multiply(*_local_field->getField());
264         if (getForcedRenormalization())
265           renormalizeTargetField(getMeasureAbsStatus());
266       }
267   }
268
269
270   /*!
271     Receives the data at time \a time in asynchronous mode. The value of the field
272     will be time-interpolated from the field values received.
273     \param time time at which the value is desired
274   */
275   void InterpKernelDEC::recvData( double time )
276   {
277     _interpolation_matrix->getAccessDEC()->setTime(time);
278     recvData() ;
279   }
280
281   /*!
282     Sends the data whether the processor is on the working side or on the lazy side.
283     It must match a recvData() call on the other side.
284   */
285   void InterpKernelDEC::sendData()
286   {
287     if (_source_group->containsMyRank())
288       {
289     
290         _interpolation_matrix->multiply(*_local_field->getField());
291         if (getForcedRenormalization())
292           renormalizeTargetField(getMeasureAbsStatus());
293     
294       }
295     else if (_target_group->containsMyRank())
296       _interpolation_matrix->transposeMultiply(*_local_field->getField());
297   }
298
299   /*!
300     Sends the data available at time \a time in asynchronous mode. 
301     \param time time at which the value is available
302     \param deltatime time interval between the value presently sent and the next one. 
303   */
304   void InterpKernelDEC::sendData( double time , double deltatime )
305   {
306     _interpolation_matrix->getAccessDEC()->setTime(time,deltatime);
307     sendData() ;
308   }
309   
310 }