5c17bb42ae7581020ec0012d399464b11ca00f7e
[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 "sci_warning.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         // check function
146         double xyz[3] = {1,1,1};
147         int numfun = 1;
148         double out;
149         double ret = int3d_f(xyz, &numfun, &out);
150         if(ret == 0)
151         {
152             Scierror(50, _("%s: Argument #%d : Variable returned by scilab argument function is incorrect.\n"), "int3d", 4);
153             DifferentialEquation::removeDifferentialEquationFunctions();
154             return types::Function::Error;
155         }
156     }
157     else if(in[3]->isString())
158     {
159         bool bOK = false;
160         types::String* pStr = in[3]->getAs<types::String>();
161         bOK = deFunctionsManager->setFFunction(pStr);
162
163         if(bOK == false)
164         {
165             char* pst = wide_string_to_UTF8(pStr->get(0));
166             Scierror(50, _("%s: Subroutine not found: %s\n"), "int3d", pst);
167             FREE(pst);
168             DifferentialEquation::removeDifferentialEquationFunctions();
169             return types::Function::Error;
170         }
171     }
172     else if(in[3]->isList())
173     {
174         types::List* pList = in[3]->getAs<types::List>();
175
176         if(pList->getSize() == 0)
177         {
178             Scierror(50, _("%s: Argument #%d : Subroutine not found in list: %s\n"), "int3d", 4, "(string empty)");
179             DifferentialEquation::removeDifferentialEquationFunctions();
180             return types::Function::Error;
181         }
182
183         if(pList->get(0)->isCallable())
184         {
185             deFunctionsManager->setFFunction(pList->get(0)->getAs<types::Callable>());
186             for(int iter = 1; iter < pList->getSize(); iter++)
187             {
188                 deFunctionsManager->setFArgs(pList->get(iter)->getAs<types::InternalType>());
189             }
190         }
191         else
192         {
193             Scierror(999, _("%s: Wrong type for input argument #%d : The first argument in the list must be a Scilab function.\n"), "int3d", 4);
194             DifferentialEquation::removeDifferentialEquationFunctions();
195             return types::Function::Error;
196         }
197     }
198     else
199     {
200         Scierror(999, _("%s: Wrong type for input argument #%d : A function expected.\n"), "int3d", 4);
201         DifferentialEquation::removeDifferentialEquationFunctions();
202         return types::Function::Error;
203     }
204
205     // nf (optional)
206     if(in.size() > 4)
207     {
208         if(in[4]->isDouble() == false)
209         {
210             Scierror(999, _("%s: Wrong type for input argument #%d : A real matrix expected.\n"), "int3d", 5);
211             DifferentialEquation::removeDifferentialEquationFunctions();
212             return types::Function::Error;
213         }
214         types::Double* pDblNf = in[4]->getAs<types::Double>();
215         if(pDblNf->isComplex())
216         {
217             Scierror(999, _("%s: Wrong type for input argument #%d : A real matrix expected.\n"), "int3d", 5);
218             DifferentialEquation::removeDifferentialEquationFunctions();
219             return types::Function::Error;
220         }
221
222         if(pDblNf->isScalar() == false)
223         {
224             Scierror(999, _("%s: Wrong size for input argument #%d : A scalar expected.\n"), "int3d", 5);
225             DifferentialEquation::removeDifferentialEquationFunctions();
226             return types::Function::Error;
227         }
228         nf = (int)pDblNf->get(0);
229
230         if(nf < 1)
231         {
232             Scierror(999, _("%s: Wrong value for input argument #%d : A positive value expected.\n"), "int3d", 5);
233             DifferentialEquation::removeDifferentialEquationFunctions();
234             return types::Function::Error;
235         }
236     }
237
238     // params (optional)
239     if(in.size() == 6)
240     {
241         if(in[5]->isDouble() == false)
242         {
243             Scierror(999, _("%s: Wrong type for input argument #%d : A real matrix expected.\n"), "int3d", 6);
244             DifferentialEquation::removeDifferentialEquationFunctions();
245             return types::Function::Error;
246         }
247
248         types::Double* pDblParams = in[5]->getAs<types::Double>();
249         if(pDblParams->isComplex())
250         {
251             Scierror(999, _("%s: Wrong type for input argument #%d : A real matrix expected.\n"), "int3d", 6);
252             DifferentialEquation::removeDifferentialEquationFunctions();
253             return types::Function::Error;
254         }
255
256         if(pDblParams->getSize() != 4)
257         {
258             Scierror(999, _("%s: Wrong size for input argument #%d : %d expected.\n"), "int3d", 6, 4);
259             DifferentialEquation::removeDifferentialEquationFunctions();
260             return types::Function::Error;
261         }
262
263         if(getWarningMode())
264         {
265             for(int i = 0; i < 4; i++)
266             {
267                 if(pDblParams->get(i) < 0)
268                 {
269                     sciprint(_("%ls: Warning : Wrong value for the element %d of argument #%d : The default value will be used.\n"),L"int3d", i+1, 6);
270                 }
271             }
272         }
273
274         minpts = pDblParams->get(0) < 0 ? minpts : (int)pDblParams->get(0);
275         maxpts = pDblParams->get(1) < 0 ? maxpts : (int)pDblParams->get(1);
276         epsabs = pDblParams->get(2) < 0.0 ? epsabs : pDblParams->get(2);
277         epsrel = pDblParams->get(3) < 0.0 ? epsrel : pDblParams->get(3);
278
279         if(pDblParams->get(2) == 0.0 && pDblParams->get(3) == 0.0)
280         {
281             if(getWarningMode())
282             {
283                 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);
284             }
285             epsabs = 0.0;
286             epsrel = 1.0e-5;
287         }
288
289         if(minpts > maxpts)
290         {
291             Scierror(999, _("%s: Wrong value for input argument #%d : minpts smaller than maxpts expected.\n"), "int3d", 6);
292             DifferentialEquation::removeDifferentialEquationFunctions();
293             return types::Function::Error;
294         }
295     }
296
297     maxpts = Max(maxpts, 43 * pDblX->getCols());
298
299 // *** Perform operation. ***
300     int one     = 1;
301     int three   = 3;
302     int maxsub  = 7 * (maxpts - 43 * pDblX->getCols()) / (8 * 43) + pDblX->getCols();
303     int size    = 4 * pDblX->getCols();
304
305     // input data : 3 rows, 4*maxsub cols
306     int dataSize = 3 * 4 * maxsub;
307     double* pdData = (double*)malloc(dataSize * sizeof(double));
308     C2F(dcopy)(&size, pDblX->get(), &one, pdData, &three);
309     C2F(dcopy)(&size, pDblY->get(), &one, pdData + 1, &three);
310     C2F(dcopy)(&size, pDblZ->get(), &one, pdData + 2, &three);
311
312     // output result
313     double* pdResult = (double*)malloc(nf * sizeof(double));
314
315     // output err
316     double* pdErr = (double*)malloc(nf * sizeof(double));
317
318     // workspace
319     int mdiv = 1; // default value, but can be changed on parallel computers. (see dcutet.f)
320     int dworkSize   = maxsub * (2 * nf + 1) + 7 * Max(8 * mdiv, pDblX->getCols()) * nf + 1;
321     int iworkSize   = maxsub + mdiv;
322     double* dwork   = (double*)malloc(dworkSize * sizeof(double));
323     int* iwork      = (int*)malloc(iworkSize * sizeof(int));
324
325     int cols = pDblX->getCols();
326     C2F(dcutet)(int3d_f, &nf, pdData, &cols, &minpts, &maxpts, &epsabs, &epsrel, &maxsub, &dworkSize, &irestar, pdResult, pdErr, &nevals, &ifail, dwork, iwork);
327
328     free(pdData);
329     free(dwork);
330     free(iwork);
331     DifferentialEquation::removeDifferentialEquationFunctions();
332
333     if(ifail)
334     {
335         if(ifail == 1)
336         {
337             if(getWarningMode())
338             {
339                 sciprint(_("%ls: Warning : maxpts was too small to obtain the required accuracy.\n"),L"int3d");
340             }
341         }
342         else if(ifail == 3)
343         {
344             Scierror(999, _("%s: The volume of one of the initially given tetrahedrons is zero.\n"), "int3d");
345             free(pdResult);
346             free(pdErr);
347             return types::Function::Error;
348         }
349         else // normaly nerver call.
350         {
351             Scierror(999, _("%s: dcutet return with error %d.\n"), "int3d", ifail);
352             free(pdResult);
353             free(pdErr);
354             return types::Function::Error;
355         }
356     }
357
358 // *** Return result in Scilab. ***
359     types::Double* pDblOut = new types::Double(nf, 1);
360     pDblOut->set(pdResult);
361     out.push_back(pDblOut);
362
363     if(_iRetCount > 1)
364     {
365         types::Double* pDblErrOut = new types::Double(nf, 1);
366         pDblErrOut->set(pdErr);
367         out.push_back(pDblErrOut);
368     }
369
370     if(_iRetCount == 3)
371     {
372         types::Double* pDblNevalsOut = new types::Double((double)nevals);
373         out.push_back(pDblNevalsOut);
374     }
375
376     free(pdResult);
377     free(pdErr);
378
379     return types::Function::OK;
380 }
381 /*--------------------------------------------------------------------------*/
382