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