source: XIOS/trunk/src/transformation/domain_algorithm_compute_connectivity.cpp @ 934

Last change on this file since 934 was 934, checked in by mhnguyen, 8 years ago

Adding new transformation: Compute_connectivity_domain

Test
+) On Curie
+) Test passes

File size: 8.1 KB
Line 
1/*!
2   \file domain_algorithm_compute_connectivity.cpp
3   \author Ha NGUYEN
4   \since 15 Jul 2016
5   \date 15 Jul 2016
6
7   \brief Algorithm for compute_connectivity on an domain.
8 */
9#include "domain_algorithm_compute_connectivity.hpp"
10#include "compute_connectivity_domain.hpp"
11#include "mesh.hpp"
12#include "domain.hpp"
13#include "grid.hpp"
14#include "grid_transformation_factory_impl.hpp"
15
16namespace xios {
17CGenericAlgorithmTransformation* CDomainAlgorithmComputeConnectivity::create(CGrid* gridDst, CGrid* gridSrc,
18                                                                     CTransformation<CDomain>* transformation,
19                                                                     int elementPositionInGrid,
20                                                                     std::map<int, int>& elementPositionInGridSrc2ScalarPosition,
21                                                                     std::map<int, int>& elementPositionInGridSrc2AxisPosition,
22                                                                     std::map<int, int>& elementPositionInGridSrc2DomainPosition,
23                                                                     std::map<int, int>& elementPositionInGridDst2ScalarPosition,
24                                                                     std::map<int, int>& elementPositionInGridDst2AxisPosition,
25                                                                     std::map<int, int>& elementPositionInGridDst2DomainPosition)
26{
27  std::vector<CDomain*> domainListDestP = gridDst->getDomains();
28  std::vector<CDomain*> domainListSrcP  = gridSrc->getDomains();
29
30  CComputeConnectivityDomain* compute_connectivityDomain = dynamic_cast<CComputeConnectivityDomain*> (transformation);
31  int domainDstIndex = elementPositionInGridDst2AxisPosition[elementPositionInGrid];
32  int domainSrcIndex = elementPositionInGridSrc2AxisPosition[elementPositionInGrid];
33
34  return (new CDomainAlgorithmComputeConnectivity(domainListDestP[domainDstIndex], domainListSrcP[domainSrcIndex], compute_connectivityDomain));
35}
36
37bool CDomainAlgorithmComputeConnectivity::registerTrans()
38{
39  CGridTransformationFactory<CDomain>::registerTransformation(TRANS_COMPUTE_CONNECTIVITY_DOMAIN, create);
40}
41
42CDomainAlgorithmComputeConnectivity::CDomainAlgorithmComputeConnectivity(CDomain* domainDestination, CDomain* domainSource,
43                                                                         CComputeConnectivityDomain* compute_connectivityDomain)
44: CDomainAlgorithmTransformation(domainDestination, domainSource)
45{
46  this->type_ = (ELEMENT_NO_MODIFICATION_WITHOUT_DATA);
47  compute_connectivityDomain->checkValid(domainDestination);
48  int& nbNeighborMax = compute_connectivityDomain->n_neighbor_max;
49  CArray<int,1>& nbNeighbor = compute_connectivityDomain->n_neighbor;
50  CArray<int,2>& localNeighbors = compute_connectivityDomain->local_neighbor;
51  switch (compute_connectivityDomain->type)
52  {
53    case CComputeConnectivityDomain::type_attr::node :
54      computeNodeConnectivity(domainDestination,
55                              nbNeighborMax,
56                              nbNeighbor,
57                              localNeighbors);
58      break;
59    case CComputeConnectivityDomain::type_attr::edge :
60      computeEdgeConnectivity(domainDestination,
61                              nbNeighborMax,
62                              nbNeighbor,
63                              localNeighbors);
64      break;
65    default:
66      break;
67  }
68}
69
70void CDomainAlgorithmComputeConnectivity::computeNodeConnectivity(CDomain* domain,
71                                                                  int& nbConnectivityMax,
72                                                                  CArray<int,1>& nbConnectivity,
73                                                                  CArray<int,2>& localConnectivity)
74{
75  CMesh mesh;
76  CArray<double,1>& lon = domain->lonvalue_1d;
77  CArray<double,1>& lat = domain->latvalue_1d;
78  CArray<double,2>& bounds_lon = domain->bounds_lon_1d;
79  CArray<double,2>& bounds_lat = domain->bounds_lat_1d;
80  mesh.createMesh(lon, lat, bounds_lon, bounds_lat);
81  CArray<int, 2>& face_nodes = mesh.face_nodes;
82  int nvertex = mesh.nvertex;
83  int nFaces  = mesh.nbFaces;
84  int nNodes  = mesh.nbNodes;
85  std::vector<std::vector<size_t> > nodeFaceConnectivity(nNodes, std::vector<size_t>());
86  for (int nf = 0; nf < nFaces; ++nf)
87  {
88    for (int nv = 0; nv < nvertex; ++nv)
89    {
90      nodeFaceConnectivity[face_nodes(nv,nf)].push_back(nf);
91    }
92  }
93
94  if (nFaces != nbConnectivity.numElements()) nbConnectivity.resize(nFaces);
95  nbConnectivityMax = 1;
96  for (int nf = 0; nf < nFaces; ++nf)
97  {
98    std::set<size_t> neighFaces;
99    for (int nv = 0; nv < nvertex; ++nv)
100    {
101      std::vector<size_t>& tmpFaces = nodeFaceConnectivity[face_nodes(nv,nf)];
102      for (int nFn = 0; nFn < tmpFaces.size(); ++nFn)
103      {
104        neighFaces.insert(tmpFaces[nFn]);
105      }
106    }
107    if (nbConnectivityMax < (neighFaces.size()-1)) nbConnectivityMax = neighFaces.size()-1;
108  }
109  if ((nbConnectivityMax != localConnectivity.shape()[0]) || (nFaces != localConnectivity.shape()[1]))
110    localConnectivity.resize(nbConnectivityMax, nFaces);
111
112  for (int nf = 0; nf < nFaces; ++nf)
113  {
114    std::set<size_t> neighFaces;
115    for (int nv = 0; nv < nvertex; ++nv)
116    {
117      std::vector<size_t>& tmpFaces = nodeFaceConnectivity[face_nodes(nv,nf)];
118      for (int nFn = 0; nFn < tmpFaces.size(); ++nFn)
119      {
120        neighFaces.insert(tmpFaces[nFn]);
121      }
122    }
123
124    neighFaces.erase(nf);
125    nbConnectivity(nf) = neighFaces.size();
126    std::set<size_t>::iterator it = neighFaces.begin(), ite = neighFaces.end();
127    for (int idx = 0; it != ite; ++it, ++idx)
128    {
129      localConnectivity(idx, nf) = *it;
130    }
131  }
132}
133
134void CDomainAlgorithmComputeConnectivity::computeEdgeConnectivity(CDomain* domain,
135                                                                  int& nbConnectivityMax,
136                                                                  CArray<int,1>& nbConnectivity,
137                                                                  CArray<int,2>& localConnectivity)
138{
139  CMesh mesh;
140  CArray<double,1>& lon = domain->lonvalue_1d;
141  CArray<double,1>& lat = domain->latvalue_1d;
142  CArray<double,2>& bounds_lon = domain->bounds_lon_1d;
143  CArray<double,2>& bounds_lat = domain->bounds_lat_1d;
144  mesh.createMesh(lon, lat, bounds_lon, bounds_lat);
145
146  CArray<int, 2>& face_edges = mesh.face_edges;
147  int nvertex = mesh.nvertex;
148  int nFaces  = mesh.nbFaces;
149  int nEdges  = mesh.nbEdges;
150  std::vector<std::vector<size_t> > edgeFaceConnectivity(nEdges, std::vector<size_t>());
151  for (int nf = 0; nf < nFaces; ++nf)
152  {
153    for (int nv = 0; nv < nvertex; ++nv)
154    {
155      edgeFaceConnectivity[face_edges(nv,nf)].push_back(nf);
156    }
157  }
158
159  if (nFaces != nbConnectivity.numElements()) nbConnectivity.resize(nFaces);
160  nbConnectivityMax = 1;
161  for (int nf = 0; nf < nFaces; ++nf)
162  {
163    std::set<size_t> neighFaces;
164    for (int nv = 0; nv < nvertex; ++nv)
165    {
166      std::vector<size_t>& tmpFaces = edgeFaceConnectivity[face_edges(nv,nf)];
167      for (int nFn = 0; nFn < tmpFaces.size(); ++nFn)
168      {
169        neighFaces.insert(tmpFaces[nFn]);
170      }
171    }
172    if (nbConnectivityMax < (neighFaces.size()-1)) nbConnectivityMax = neighFaces.size() - 1;
173  }
174  if ((nbConnectivityMax != localConnectivity.shape()[0]) || (nFaces != localConnectivity.shape()[1]))
175    localConnectivity.resize(nbConnectivityMax, nFaces);
176
177  for (int nf = 0; nf < nFaces; ++nf)
178  {
179    std::set<size_t> neighFaces;
180    for (int nv = 0; nv < nvertex; ++nv)
181    {
182      std::vector<size_t>& tmpFaces = edgeFaceConnectivity[face_edges(nv,nf)];
183      for (int nFn = 0; nFn < tmpFaces.size(); ++nFn)
184      {
185        neighFaces.insert(tmpFaces[nFn]);
186      }
187    }
188
189    neighFaces.erase(nf);
190    nbConnectivity(nf) = neighFaces.size();
191    std::set<size_t>::iterator it = neighFaces.begin(), ite = neighFaces.end();
192    for (int idx = 0; it != ite; ++it, ++idx)
193    {
194      localConnectivity(idx, nf) = *it;
195    }
196  }
197
198}
199
200/*!
201  Compute the index mapping between domain on grid source and one on grid destination
202*/
203void CDomainAlgorithmComputeConnectivity::computeIndexSourceMapping_(const std::vector<CArray<double,1>* >& dataAuxInputs)
204{
205
206}
207
208}
Note: See TracBrowser for help on using the repository browser.