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