[interpolation] mesh2d gateway introduced
[scilab.git] / scilab / modules / interpolation / sci_gateway / cpp / sci_mesh2di.cpp
1 /*
2  *  Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
3  *  Copyright (C) 2018 - St├ęphane MOTTELET
4  *
5  * This file is hereby licensed under the terms of the GNU GPL v2.0,
6  * For more information, see the COPYING file which you should have received
7  * along with this program.
8  *
9  */
10 /*--------------------------------------------------------------------------*/
11 #include "interpolation_gw.hxx"
12 #include "function.hxx"
13 #include "double.hxx"
14 #include "sparse.hxx"
15 #include <numeric>
16
17
18 extern "C"
19 {
20 #include "localization.h"
21 #include "Scierror.h"
22 #include "interpolation_functions.h"
23 }
24 /*--------------------------------------------------------------------------*/
25
26 types::Function::ReturnValue sci_mesh2di(types::typed_list &in, int _iRetCount, types::typed_list &out)
27 {
28     if (in.size() != 3)
29     {
30         Scierror(77, _("%s: Wrong number of input argument(s): %d expected.\n"), "mesh2di", 3);
31         return types::Function::Error;
32     }
33
34     if (_iRetCount > 1)
35     {
36         Scierror(78, _("%s: Wrong number of output argument(s): %d expected.\n"), "mesh2di", 1);
37         return types::Function::Error;
38     }
39
40     for (int i = 0; i < in.size(); i++)
41     {
42         if (!in[i]->isDouble() | in[i]->getAs<types::Double>()->isComplex())
43         {
44             Scierror(999, _("%s: Wrong type for argument #%d: Real matrix expected.\n"), "mesh2di", i);
45             return types::Function::Error;
46         }
47     }
48
49     types::Double *pDblx = in[0]->getAs<types::Double>();
50     types::Double *pDbly = in[1]->getAs<types::Double>();
51
52     if (pDblx->getSize() != pDblx->getSize())
53     {
54         Scierror(999, _("%s: Arguments #%d and #%d must have the same sizes.\n"), "mesh2di", 1,2);
55         return types::Function::Error;
56     }
57
58     int iNbNodes = pDblx->getSize();
59
60     if (iNbNodes < 3)
61     {
62         Scierror(999, _("%s: Wrong size for input argument #%d : At least %d elements expected.\n"), "mesh2di", 1, 3);
63         return types::Function::Error;
64     }
65     
66     double* pdblx = pDblx->get();
67     double* pdbly = pDbly->get();
68     std::vector<int> bdy;
69
70     // verify boundary specification
71     types::Double* pDblBdy = in[2]->getAs<types::Double>();
72     int iStart = -1;
73
74     for (int i=0; i < pDblBdy->getSize(); i++)
75     {
76         double dblCurrNode = pDblBdy->get(i);
77         if (dblCurrNode != std::round(dblCurrNode) || (int)dblCurrNode < 1 || (int)dblCurrNode > iNbNodes)
78         {
79             Scierror(999, _("%s: Wrong value for input argument #%d: incorrect element #%d\n"), "mesh2di", 3, i+1);
80             return types::Function::Error;
81         }
82         // if (iStart < 0) next node is the first of a connected component
83         // else if iCurrNode == iStart, close the current component (set iStart = -1)
84         iStart = iStart < 0 ? (int)dblCurrNode : ((int)dblCurrNode == iStart ? -1 : iStart);
85         bdy.push_back((int)dblCurrNode);
86     }
87     // iStart should be equal to -1 (all components closed)
88     if (iStart != -1)
89     {
90         Scierror(999, _("%s: last connected component of boundary is not closed\n"), "mesh2di");
91         return types::Function::Error;
92     }
93
94     // allocate working arrays
95     int* piTri = NULL;
96     int iNbTri = 0;
97     int iErr = 0;
98     int in4 = 4*iNbNodes-4;
99     int in6 = 6*(2*iNbNodes-2);
100     int* piWork = new int[in4];
101     piTri = new int[3*in4];
102     int iNbBdyNodes = bdy.size();
103     
104     // mesh2b wants xy=[x(:)';y(:)'] + a clone of it
105     types::Double* pDblxy = new types::Double(2,iNbNodes);
106     for (int i=0; i < iNbNodes; i++)
107     {
108          pDblxy->set(0,i,pDblx->get(i));
109          pDblxy->set(1,i,pDbly->get(i));
110     }
111     types::Double* pDblxyClone = pDblxy->clone();
112
113     C2F(mesh2b)(&iNbNodes, &in6, &in4, &iNbBdyNodes, pDblxy->get(), pDblxyClone->get(),
114                  piTri, piWork, &bdy.front(), &iNbTri, &iErr);
115
116     pDblxy->killMe();
117     pDblxyClone->killMe();
118     delete[] piWork;
119
120     if (iErr)
121     {
122         std::string errStr[] = {
123         _("%s: some points are identical.\n"),
124         _("%s: some points are identical.\n"),
125         _("%s: all points are aligned.\n"),
126         _("%s: wrong boundary array.\n"),
127         _("%s: crossed boundary.\n"),
128         _("%s: wrong orientation of the boundary.\n"),
129         _("%s: size limitation.\n"),
130         _("%s: crossed boundary.\n"),
131         _("%s: an interior point is on the boundary.\n"),
132         _("%s: an interior point is too close to the boundary.\n"),
133         _("%s: some points are identical or size limitation.\n")};
134
135         Scierror(999, errStr[iErr-1].c_str(), "mesh2di");
136
137         delete[] piTri;
138
139         return types::Function::Error;
140     }
141     
142     // fix empty triangles + eventual wrong orientation (should be counter clock-wise)
143     types::Double* pDblTrianglesNodes = new types::Double(3,iNbTri);
144     double* pdblTri = pDblTrianglesNodes->get();
145     int iTriLin = 0;
146     int iTriLinOrig = 0;
147     int i0, i1, i2;
148
149     for (int i = 0; i < iNbTri; i++)
150     {
151         i0 = piTri[iTriLinOrig++];
152         i1 = piTri[iTriLinOrig++];
153         i2 = piTri[iTriLinOrig++];
154         if (i0>0) // triangle is not empty
155         {
156             pdblTri[iTriLin++] = (double)i0;
157             // compute determinant of (i0->i1) and (i0->i2) vectors
158             if ((pdblx[i1-1]-pdblx[i0-1])*(pdbly[i2-1]-pdbly[i0-1])-
159                 (pdblx[i2-1]-pdblx[i0-1])*(pdbly[i1-1]-pdbly[i0-1]) > 0.0)
160             {
161                 pdblTri[iTriLin++] = (double)i1;
162                 pdblTri[iTriLin++] = (double)i2;
163             }
164             else
165             {
166                 pdblTri[iTriLin++] = (double)i2;
167                 pdblTri[iTriLin++] = (double)i1;
168             }
169         }
170     }
171     delete[] piTri;
172     pDblTrianglesNodes->resize(3,iTriLin/3);
173
174     out.push_back(pDblTrianglesNodes);
175
176     return types::Function::OK;
177 }
178