Merge remote-tracking branch 'origin/master' into windows
[scilab.git] / scilab / modules / differential_equations / sci_gateway / cpp / sci_int3d.cpp
1 /*
2 * Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
3 * Copyright (C) 2011 - DIGITEO - Cedric DELAMARRE
4 *
5  * Copyright (C) 2012 - 2016 - Scilab Enterprises
6  *
7  * This file is hereby licensed under the terms of the GNU GPL v2.0,
8  * pursuant to article 5.3.4 of the CeCILL v.2.1.
9  * This file was originally licensed under the terms of the CeCILL v2.1,
10  * and continues to be available under such terms.
11  * For more information, see the COPYING file which you should have received
12  * along with this program.
13 *
14 */
15 /*--------------------------------------------------------------------------*/
16
17 #include "differential_equations_gw.hxx"
18 #include "function.hxx"
19 #include "double.hxx"
20 #include "string.hxx"
21 #include "list.hxx"
22 #include "callable.hxx"
23 #include "differentialequationfunctions.hxx"
24
25 extern "C"
26 {
27 #include "sci_malloc.h"
28 #include "elem_common.h"
29 #include "localization.h"
30 #include "Scierror.h"
31 #include "scifunctions.h"
32 #include "configvariable_interface.h"
33 #include "sciprint.h"
34 }
35
36 /*--------------------------------------------------------------------------*/
37 types::Function::ReturnValue sci_int3d(types::typed_list &in, int _iRetCount, types::typed_list &out)
38 {
39     //input
40     types::Double* pDblX = NULL;
41     types::Double* pDblY = NULL;
42     types::Double* pDblZ = NULL;
43
44     int nf          = 1;
45     int minpts      = 0;
46     int maxpts      = 1000;
47     double epsabs   = 0.0;
48     double epsrel   = 1.0e-5;
49     int irestar     = 0;
50
51     // output
52     int ifail = 0;
53     int nevals = 0;
54
55     // error message catched
56     std::ostringstream os;
57     bool bCatch = false;
58
59     // *** check the minimal number of input args. ***
60     if (in.size() < 4 || in.size() > 6)
61     {
62         Scierror(77, _("%s: Wrong number of input argument(s): %d to %d expected.\n"), "int3d", 4, 6);
63         return types::Function::Error;
64     }
65
66     // *** check number of output args according the methode. ***
67     if (_iRetCount > 3)
68     {
69         Scierror(78, _("%s: Wrong number of output argument(s): %d expected.\n"), "int3d", 2);
70         return types::Function::Error;
71     }
72
73     // *** check type of input args and get it. ***
74     // X
75     if (in[0]->isDouble() == false)
76     {
77         Scierror(999, _("%s: Wrong type for input argument #%d: A real matrix expected.\n"), "int3d", 1);
78         return types::Function::Error;
79     }
80     pDblX = in[0]->getAs<types::Double>();
81     if (pDblX->isComplex())
82     {
83         Scierror(999, _("%s: Wrong type for input argument #%d: A real matrix expected.\n"), "int3d", 1);
84         return types::Function::Error;
85     }
86
87     if (pDblX->getRows() != 4)
88     {
89         Scierror(999, _("%s: Wrong size for input argument #%d: A 4 by N matrix expected.\n"), "int3d", 1);
90         return types::Function::Error;
91     }
92
93     // Y
94     if (in[1]->isDouble() == false)
95     {
96         Scierror(999, _("%s: Wrong type for input argument #%d: A real matrix expected.\n"), "int3d", 2);
97         return types::Function::Error;
98     }
99     pDblY = in[1]->getAs<types::Double>();
100     if (pDblY->isComplex())
101     {
102         Scierror(999, _("%s: Wrong type for input argument #%d: A real matrix expected.\n"), "int3d", 2);
103         return types::Function::Error;
104     }
105
106     if (pDblY->getRows() != 4)
107     {
108         Scierror(999, _("%s: Wrong size for input argument #%d: A 4 by N matrix expected.\n"), "int3d", 2);
109         return types::Function::Error;
110     }
111
112     if (pDblY->getCols() != pDblX->getCols())
113     {
114         Scierror(999, _("%s: Wrong size for input argument #%d: Same size of input argument %d expected.\n"), "int3d", 2, 1);
115         return types::Function::Error;
116     }
117
118     // Z
119     if (in[2]->isDouble() == false)
120     {
121         Scierror(999, _("%s: Wrong type for input argument #%d: A real matrix expected.\n"), "int3d", 3);
122         return types::Function::Error;
123     }
124     pDblZ = in[2]->getAs<types::Double>();
125     if (pDblZ->isComplex())
126     {
127         Scierror(999, _("%s: Wrong type for input argument #%d: A real matrix expected.\n"), "int3d", 3);
128         return types::Function::Error;
129     }
130
131     if (pDblZ->getRows() != 4)
132     {
133         Scierror(999, _("%s: Wrong size for input argument #%d: A 4 by N matrix expected.\n"), "int3d", 4);
134         return types::Function::Error;
135     }
136
137     if (pDblZ->getCols() != pDblX->getCols())
138     {
139         Scierror(999, _("%s: Wrong size for input argument #%d: Same size of input argument %d expected.\n"), "int3d", 3, 1);
140         return types::Function::Error;
141     }
142
143     // function
144     DifferentialEquationFunctions deFunctionsManager("int3d");
145     DifferentialEquation::addDifferentialEquationFunctions(&deFunctionsManager);
146
147     if (in[3]->isCallable())
148     {
149         types::Callable* pCall = in[3]->getAs<types::Callable>();
150         deFunctionsManager.setFFunction(pCall);
151     }
152     else if (in[3]->isString())
153     {
154         bool bOK = false;
155         types::String* pStr = in[3]->getAs<types::String>();
156         bOK = deFunctionsManager.setFFunction(pStr);
157
158         if (bOK == false)
159         {
160             const char* pst = pStr->get(0);
161             Scierror(50, _("%s: Subroutine not found: %s\n"), "int3d", pst);
162             DifferentialEquation::removeDifferentialEquationFunctions();
163             return types::Function::Error;
164         }
165     }
166     else if (in[3]->isList())
167     {
168         types::List* pList = in[3]->getAs<types::List>();
169
170         if (pList->getSize() == 0)
171         {
172             Scierror(50, _("%s: Argument #%d: Subroutine not found in list: %s\n"), "int3d", 4, "(string empty)");
173             DifferentialEquation::removeDifferentialEquationFunctions();
174             return types::Function::Error;
175         }
176
177         if (pList->get(0)->isCallable())
178         {
179             deFunctionsManager.setFFunction(pList->get(0)->getAs<types::Callable>());
180             for (int iter = 1; iter < pList->getSize(); iter++)
181             {
182                 deFunctionsManager.setFArgs(pList->get(iter)->getAs<types::InternalType>());
183             }
184         }
185         else
186         {
187             Scierror(999, _("%s: Wrong type for input argument #%d: The first argument in the list must be a Scilab function.\n"), "int3d", 4);
188             DifferentialEquation::removeDifferentialEquationFunctions();
189             return types::Function::Error;
190         }
191     }
192     else
193     {
194         Scierror(999, _("%s: Wrong type for input argument #%d: A function expected.\n"), "int3d", 4);
195         DifferentialEquation::removeDifferentialEquationFunctions();
196         return types::Function::Error;
197     }
198
199     // nf (optional)
200     if (in.size() > 4)
201     {
202         if (in[4]->isDouble() == false)
203         {
204             Scierror(999, _("%s: Wrong type for input argument #%d: A real matrix expected.\n"), "int3d", 5);
205             DifferentialEquation::removeDifferentialEquationFunctions();
206             return types::Function::Error;
207         }
208         types::Double* pDblNf = in[4]->getAs<types::Double>();
209         if (pDblNf->isComplex())
210         {
211             Scierror(999, _("%s: Wrong type for input argument #%d: A real matrix expected.\n"), "int3d", 5);
212             DifferentialEquation::removeDifferentialEquationFunctions();
213             return types::Function::Error;
214         }
215
216         if (pDblNf->isScalar() == false)
217         {
218             Scierror(999, _("%s: Wrong size for input argument #%d: A scalar expected.\n"), "int3d", 5);
219             DifferentialEquation::removeDifferentialEquationFunctions();
220             return types::Function::Error;
221         }
222         nf = (int)pDblNf->get(0);
223
224         if (nf < 1)
225         {
226             Scierror(999, _("%s: Wrong value for input argument #%d: A positive value expected.\n"), "int3d", 5);
227             DifferentialEquation::removeDifferentialEquationFunctions();
228             return types::Function::Error;
229         }
230     }
231
232     // params (optional)
233     if (in.size() == 6)
234     {
235         if (in[5]->isDouble() == false)
236         {
237             Scierror(999, _("%s: Wrong type for input argument #%d: A real matrix expected.\n"), "int3d", 6);
238             DifferentialEquation::removeDifferentialEquationFunctions();
239             return types::Function::Error;
240         }
241
242         types::Double* pDblParams = in[5]->getAs<types::Double>();
243         if (pDblParams->isComplex())
244         {
245             Scierror(999, _("%s: Wrong type for input argument #%d: A real matrix expected.\n"), "int3d", 6);
246             DifferentialEquation::removeDifferentialEquationFunctions();
247             return types::Function::Error;
248         }
249
250         if (pDblParams->getSize() != 4)
251         {
252             Scierror(999, _("%s: Wrong size for input argument #%d: %d expected.\n"), "int3d", 6, 4);
253             DifferentialEquation::removeDifferentialEquationFunctions();
254             return types::Function::Error;
255         }
256
257         if (getWarningMode())
258         {
259             for (int i = 0; i < 4; i++)
260             {
261                 if (pDblParams->get(i) < 0)
262                 {
263                     sciprint(_("%s: Warning: Wrong value for the element %d of argument #%d: The default value will be used.\n"), "int3d", i + 1, 6);
264                 }
265             }
266         }
267
268         minpts = pDblParams->get(0) < 0 ? minpts : (int)pDblParams->get(0);
269         maxpts = pDblParams->get(1) < 0 ? maxpts : (int)pDblParams->get(1);
270         epsabs = pDblParams->get(2) < 0.0 ? epsabs : pDblParams->get(2);
271         epsrel = pDblParams->get(3) < 0.0 ? epsrel : pDblParams->get(3);
272
273         if (pDblParams->get(2) == 0.0 && pDblParams->get(3) == 0.0)
274         {
275             if (getWarningMode())
276             {
277                 sciprint(_("%s: Warning: Wrong value for the element %d and %d of argument #%d: The default value will be used.\n"), "int3d", 3, 4, 6);
278             }
279             epsabs = 0.0;
280             epsrel = 1.0e-5;
281         }
282
283         if (minpts > maxpts)
284         {
285             Scierror(999, _("%s: Wrong value for input argument #%d: minpts smaller than maxpts expected.\n"), "int3d", 6);
286             DifferentialEquation::removeDifferentialEquationFunctions();
287             return types::Function::Error;
288         }
289     }
290
291     maxpts = std::max(maxpts, 43 * pDblX->getCols());
292
293     // *** Perform operation. ***
294     int one     = 1;
295     int three   = 3;
296     int maxsub  = 7 * (maxpts - 43 * pDblX->getCols()) / (8 * 43) + pDblX->getCols();
297     int size    = 4 * pDblX->getCols();
298
299     // input data: 3 rows, 4*maxsub cols
300     int dataSize = 3 * 4 * maxsub;
301     double* pdData = (double*)MALLOC(dataSize * sizeof(double));
302     C2F(dcopy)(&size, pDblX->get(), &one, pdData, &three);
303     C2F(dcopy)(&size, pDblY->get(), &one, pdData + 1, &three);
304     C2F(dcopy)(&size, pDblZ->get(), &one, pdData + 2, &three);
305
306     // output result
307     double* pdResult = (double*)MALLOC(nf * sizeof(double));
308
309     // output err
310     double* pdErr = (double*)MALLOC(nf * sizeof(double));
311
312     // workspace
313     int mdiv = 1; // default value, but can be changed on parallel computers. (see dcutet.f)
314     int dworkSize   = maxsub * (2 * nf + 1) + 7 * std::max(8 * mdiv, pDblX->getCols()) * nf + 1;
315     int iworkSize   = maxsub + mdiv;
316     double* dwork   = (double*)MALLOC(dworkSize * sizeof(double));
317     int* iwork      = (int*)MALLOC(iworkSize * sizeof(int));
318
319     int cols = pDblX->getCols();
320     try
321     {
322         C2F(dcutet)(int3d_f, &nf, pdData, &cols, &minpts, &maxpts, &epsabs, &epsrel, &maxsub, &dworkSize, &irestar, pdResult, pdErr, &nevals, &ifail, dwork, iwork);
323     }
324     catch (ast::InternalError &ie)
325     {
326         os << ie.GetErrorMessage();
327         bCatch = true;
328     }
329
330     FREE(pdData);
331     FREE(dwork);
332     FREE(iwork);
333     DifferentialEquation::removeDifferentialEquationFunctions();
334
335     if (bCatch)
336     {
337         char szError[bsiz];
338         os_sprintf(szError, bsiz, _("%s: An error occured in '%s' subroutine.\n"), "int3d", "dcutet");
339         os << szError;
340         throw ast::InternalError(os.str());
341     }
342
343     if (ifail)
344     {
345         if (ifail == 1)
346         {
347             if (getWarningMode())
348             {
349                 sciprint(_("%s: Warning: maxpts was too small to obtain the required accuracy.\n"), "int3d");
350             }
351         }
352         else if (ifail == 3)
353         {
354             Scierror(999, _("%s: The volume of one of the initially given tetrahedrons is zero.\n"), "int3d");
355             FREE(pdResult);
356             FREE(pdErr);
357             return types::Function::Error;
358         }
359         else // normaly nerver call.
360         {
361             Scierror(999, _("%s: dcutet return with error %d.\n"), "int3d", ifail);
362             FREE(pdResult);
363             FREE(pdErr);
364             return types::Function::Error;
365         }
366     }
367
368     // *** Return result in Scilab. ***
369     types::Double* pDblOut = new types::Double(nf, 1);
370     pDblOut->set(pdResult);
371     out.push_back(pDblOut);
372
373     if (_iRetCount > 1)
374     {
375         types::Double* pDblErrOut = new types::Double(nf, 1);
376         pDblErrOut->set(pdErr);
377         out.push_back(pDblErrOut);
378     }
379
380     if (_iRetCount == 3)
381     {
382         types::Double* pDblNevalsOut = new types::Double((double)nevals);
383         out.push_back(pDblNevalsOut);
384     }
385
386     FREE(pdResult);
387     FREE(pdErr);
388
389     return types::Function::OK;
390 }
391 /*--------------------------------------------------------------------------*/
392