]> SALOME platform Git repositories - tools/medcoupling.git/blob - src/ParaMEDMEM/InterpKernelDEC.cxx
Salome HOME
MERGE stage 1: keep doc/dev and src/MEDCalc/doc
[tools/medcoupling.git] / src / ParaMEDMEM / InterpKernelDEC.cxx
1 // Copyright (C) 2007-2015  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 ParaMEDMEM
34 {  
35
36   /*!
37     \anchor InterpKernelDEC-det
38     \class InterpKernelDEC
39
40     \section dec-over Overview
41
42     The InterpKernelDEC enables the \ref InterpKerRemapGlobal "remapping" of fields between two parallel codes.
43     This remapping is based on the computation of intersection volumes between elements from code A
44     and elements from code B. The computation is possible for 3D meshes, 2D meshes, and 3D-surface
45     meshes. Dimensions must be similar for code A and code B (for instance, though it could be
46     desirable, it is not yet possible to couple 3D surfaces with 2D surfaces).
47
48     In the present version, only fields lying on elements are considered.
49
50     \image html NonCoincident_small.png "Example showing the transfer from a field based on a
51     quadrangular mesh to a triangular mesh. In a P0-P0 interpolation, to obtain the value on a triangle,
52     the values on quadrangles are weighted by their intersection area and summed."
53
54     \image latex NonCoincident_small.eps "Example showing the transfer from a field based on a quadrangular
55      mesh to a triangular mesh. In a P0-P0 interpolation, to obtain the value on a triangle, the values
56      on quadrangles are weighted by their intersection area and summed."
57
58     A typical use of InterpKernelDEC encompasses two distinct phases :
59     - A setup phase during which the intersection volumes are computed and the communication structures are
60     setup. This corresponds to calling the InterpKernelDEC::synchronize() method.
61     - A use phase during which the remappings are actually performed. This corresponds to the calls to
62     sendData() and recvData() which actually trigger the data exchange. The data exchange are synchronous
63     in the current version of the library so that recvData() and sendData() calls must be synchronized
64     on code A and code B processor groups.
65
66     The following code excerpt illutrates a typical use of the InterpKernelDEC class.
67
68     \code
69     ...
70     InterpKernelDEC dec(groupA, groupB);
71     dec.attachLocalField(field);
72     dec.synchronize();
73     if (groupA.containsMyRank())
74     dec.recvData();
75     else if (groupB.containsMyRank())
76     dec.sendData();
77     ...
78     \endcode
79     A \ref InterpKerRemapGlobal "remapping" of the field from the source mesh to the target mesh is performed by
80     the function synchronise(), which computes the interpolation matrix.
81
82     Computing the field on the receiving side can be expressed in terms of a matrix-vector product :
83     \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
84     on the source side.
85     When remapping a 3D surface to another 3D surface, a projection phase is necessary to match elements
86     from both sides. Care must be taken when defining this projection to obtain a
87     \ref InterpKerRemapGlobal "conservative remapping".
88
89     In the P0-P0 case, this matrix is a plain rectangular matrix with coefficients equal to the
90     intersection areas between triangle and quadrangles. For instance, in the above figure, the matrix
91     is :
92
93     \f[
94     \begin{tabular}{|cccc|}
95     0.72 & 0 & 0.2 & 0 \\
96     0.46 & 0 & 0.51 & 0.03\\
97     0.42 & 0.53 & 0 & 0.05\\
98     0 & 0 & 0.92 & 0.05 \\
99     \end{tabular}
100     \f]
101
102
103
104     \section interpkerneldec_options Options
105     On top of \ref dec_options, options supported by %InterpKernelDEC objects are
106     related to the underlying Intersector class. 
107     All the options available in the intersector objects are
108     available for the %InterpKernelDEC object. The various options available for  * intersectors can
109     be reviewed in \ref InterpKerIntersectors.
110  
111     For instance :
112     \verbatim
113     InterpKernelDEC dec(source_group, target_group);
114     dec.attachLocalField(field);
115     dec.setOptions("DoRotate",false);
116     dec.setOptions("Precision",1e-12);
117     dec.synchronize();
118     \endverbatim
119
120     \warning{  Options must be set before calling the synchronize method. }
121   */
122   
123   InterpKernelDEC::InterpKernelDEC():_interpolation_matrix(0)
124   {  
125   }
126
127   /*!
128     This constructor creates an InterpKernelDEC which has \a source_group as a working side 
129     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.
130     The constructor must be called synchronously on all processors of both processor groups.
131
132     \param source_group working side ProcessorGroup
133     \param target_group lazy side ProcessorGroup
134
135   */
136   InterpKernelDEC::InterpKernelDEC(ProcessorGroup& source_group, ProcessorGroup& target_group):
137     DisjointDEC(source_group, target_group),_interpolation_matrix(0)
138   {
139
140   }
141
142   InterpKernelDEC::InterpKernelDEC(const std::set<int>& src_ids, const std::set<int>& trg_ids,
143                                    const MPI_Comm& world_comm):DisjointDEC(src_ids,trg_ids,world_comm),
144                                                                _interpolation_matrix(0)
145   {
146   }
147
148   InterpKernelDEC::~InterpKernelDEC()
149   {
150     if (_interpolation_matrix !=0)
151       delete _interpolation_matrix;
152   } 
153
154   /*! 
155     \brief Synchronization process for exchanging topologies.
156
157     This method prepares all the structures necessary for sending data from a processor group to the other. It uses the mesh underlying the fields that have been set with attachLocalField method.
158     It works in four steps :
159     -# Bounding boxes are computed for each subdomain,
160     -# The lazy side mesh parts that are likely to intersect the working side local processor are sent to the working side,
161     -# The working side calls the interpolation kernel to compute the intersection between local and imported mesh.
162     -# The lazy side is updated so that it knows the structure of the data that will be sent by
163     the working side during a \a sendData() call.
164
165   */
166   void InterpKernelDEC::synchronize()
167   {
168     if(!isInUnion())
169       return ;
170     delete _interpolation_matrix;
171     _interpolation_matrix = new InterpolationMatrix (_local_field, *_source_group,*_target_group,*this,*this); 
172
173     //setting up the communication DEC on both sides  
174     if (_source_group->containsMyRank())
175       {
176         //locate the distant meshes
177         ElementLocator locator(*_local_field, *_target_group, *_source_group);
178         //transfering option from InterpKernelDEC to ElementLocator   
179         locator.copyOptions(*this);
180         MEDCouplingPointSet* distant_mesh=0; 
181         int* distant_ids=0;
182         std::string distantMeth;
183         for (int i=0; i<_target_group->size(); i++)
184           {
185             //        int idistant_proc = (i+_source_group->myRank())%_target_group->size();
186             int idistant_proc=i;
187
188             //gathers pieces of the target meshes that can intersect the local mesh
189             locator.exchangeMesh(idistant_proc,distant_mesh,distant_ids);
190             if (distant_mesh !=0)
191               {
192                 locator.exchangeMethod(_method,idistant_proc,distantMeth);
193                 //adds the contribution of the distant mesh on the local one
194                 int idistant_proc_in_union=_union_group->translateRank(_target_group,idistant_proc);
195                 //std::cout <<"add contribution from proc "<<idistant_proc_in_union<<" to proc "<<_union_group->myRank()<<std::endl;
196                 _interpolation_matrix->addContribution(*distant_mesh,idistant_proc_in_union,distant_ids,_method,distantMeth);
197                 distant_mesh->decrRef();
198                 delete [] distant_ids;
199                 distant_mesh=0;
200                 distant_ids=0;
201               }
202           }
203        _interpolation_matrix->finishContributionW(locator);
204       }
205
206     if (_target_group->containsMyRank())
207       {
208         ElementLocator locator(*_local_field, *_source_group, *_target_group);
209         //transfering option from InterpKernelDEC to ElementLocator
210         locator.copyOptions(*this);
211         MEDCouplingPointSet* distant_mesh=0;
212         int* distant_ids=0;
213         for (int i=0; i<_source_group->size(); i++)
214           {
215             //        int idistant_proc = (i+_target_group->myRank())%_source_group->size();
216             int  idistant_proc=i;
217             //gathers pieces of the target meshes that can intersect the local mesh
218             locator.exchangeMesh(idistant_proc,distant_mesh,distant_ids);
219             //std::cout << " Data sent from "<<_union_group->myRank()<<" to source proc "<< idistant_proc<<std::endl;
220             if (distant_mesh!=0)
221               {
222                 std::string distantMeth;
223                 locator.exchangeMethod(_method,idistant_proc,distantMeth);
224                 distant_mesh->decrRef();
225                 delete [] distant_ids;
226                 distant_mesh=0;
227                 distant_ids=0;
228               }
229           }
230         _interpolation_matrix->finishContributionL(locator);
231       }
232     _interpolation_matrix->prepare();
233   }
234
235
236   /*!
237     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.
238   */
239   void InterpKernelDEC::recvData()
240   {
241     if (_source_group->containsMyRank())
242       _interpolation_matrix->transposeMultiply(*_local_field->getField());
243     else if (_target_group->containsMyRank())
244       {
245         _interpolation_matrix->multiply(*_local_field->getField());
246         if (getForcedRenormalization())
247           renormalizeTargetField(getMeasureAbsStatus());
248       }
249   }
250
251
252   /*!
253     Receives the data at time \a time in asynchronous mode. The value of the field
254     will be time-interpolated from the field values received.
255     \param time time at which the value is desired
256   */
257   void InterpKernelDEC::recvData( double time )
258   {
259     _interpolation_matrix->getAccessDEC()->setTime(time);
260     recvData() ;
261   }
262
263   /*!
264     Sends the data whether the processor is on the working side or on the lazy side.
265     It must match a recvData() call on the other side.
266   */
267   void InterpKernelDEC::sendData()
268   {
269     if (_source_group->containsMyRank())
270       {
271     
272         _interpolation_matrix->multiply(*_local_field->getField());
273         if (getForcedRenormalization())
274           renormalizeTargetField(getMeasureAbsStatus());
275     
276       }
277     else if (_target_group->containsMyRank())
278       _interpolation_matrix->transposeMultiply(*_local_field->getField());
279   }
280
281   /*!
282     Sends the data available at time \a time in asynchronous mode. 
283     \param time time at which the value is available
284     \param deltatime time interval between the value presently sent and the next one. 
285   */
286   void InterpKernelDEC::sendData( double time , double deltatime )
287   {
288     _interpolation_matrix->getAccessDEC()->setTime(time,deltatime);
289     sendData() ;
290   }
291   
292 }