0022786: 64bit issue in AdvApp2Var
[occt.git] / src / AdvApp2Var / AdvApp2Var_ApproxF2var.cxx
1 // Copyright (c) 1999-2012 OPEN CASCADE SAS
2 //
3 // The content of this file is subject to the Open CASCADE Technology Public
4 // License Version 6.5 (the "License"). You may not use the content of this file
5 // except in compliance with the License. Please obtain a copy of the License
6 // at http://www.opencascade.org and read it completely before using this file.
7 //
8 // The Initial Developer of the Original Code is Open CASCADE S.A.S., having its
9 // main offices at: 1, place des Freres Montgolfier, 78280 Guyancourt, France.
10 //
11 // The Original Code and all software distributed under the License is
12 // distributed on an "AS IS" basis, without warranty of any kind, and the
13 // Initial Developer hereby disclaims all such warranties, including without
14 // limitation, any warranties of merchantability, fitness for a particular
15 // purpose or non-infringement. Please see the License for the specific terms
16 // and conditions governing the rights and limitations under the License.
17
18 // AdvApp2Var_ApproxF2var.cxx
19 #include <math.h>
20 #include <AdvApp2Var_SysBase.hxx>
21 #include <AdvApp2Var_MathBase.hxx>
22 #include <AdvApp2Var_Data_f2c.hxx>
23 #include <AdvApp2Var_Data.hxx>
24 #include <AdvApp2Var_ApproxF2var.hxx>
25
26
27 static
28 int mmjacpt_(const integer *ndimen,
29              const integer *ncoefu, 
30              const integer *ncoefv, 
31              const integer *iordru, 
32              const integer *iordrv, 
33              const doublereal *ptclgd, 
34              doublereal *ptcaux, 
35              doublereal *ptccan);
36
37
38
39 static
40 int mma2ce2_(integer *numdec, 
41              integer *ndimen, 
42              integer *nbsesp, 
43              integer *ndimse, 
44              integer *ndminu, 
45              integer *ndminv, 
46              integer *ndguli, 
47              integer *ndgvli, 
48              integer *ndjacu, 
49              integer *ndjacv, 
50              integer *iordru, 
51              integer *iordrv, 
52              integer *nbpntu, 
53              integer *nbpntv, 
54              doublereal *epsapr, 
55              doublereal *sosotb, 
56              doublereal *disotb, 
57              doublereal *soditb, 
58              doublereal *diditb, 
59              doublereal *gssutb, 
60              doublereal *gssvtb, 
61              doublereal *xmaxju, 
62              doublereal *xmaxjv, 
63              doublereal *vecerr, 
64              doublereal *chpair, 
65              doublereal *chimpr, 
66              doublereal *patjac, 
67              doublereal *errmax, 
68              doublereal *errmoy, 
69              integer *ndegpu, 
70              integer *ndegpv, 
71              integer *itydec, 
72              integer *iercod);
73
74 static
75 int mma2cfu_(integer *ndujac, 
76              integer *nbpntu, 
77              integer *nbpntv, 
78              doublereal *sosotb, 
79              doublereal *disotb, 
80              doublereal *soditb, 
81              doublereal *diditb, 
82              doublereal *gssutb, 
83              doublereal *chpair, 
84              doublereal *chimpr);
85
86 static
87 int mma2cfv_(integer *ndvjac, 
88              integer *mindgu,
89              integer *maxdgu, 
90              integer *nbpntv, 
91              doublereal *gssvtb, 
92              doublereal *chpair, 
93              doublereal *chimpr, 
94              doublereal *patjac);
95
96 static
97 int mma2er1_(integer *ndjacu, 
98              integer *ndjacv, 
99              integer *ndimen, 
100              integer *mindgu, 
101              integer *maxdgu, 
102              integer *mindgv, 
103              integer *maxdgv, 
104              integer *iordru, 
105              integer *iordrv, 
106              doublereal *xmaxju, 
107              doublereal *xmaxjv, 
108              doublereal *patjac, 
109              doublereal *vecerr, 
110              doublereal *erreur);
111
112 static
113 int mma2er2_(integer *ndjacu, 
114              integer *ndjacv,
115              integer *ndimen, 
116              integer *mindgu, 
117              integer *maxdgu, 
118              integer *mindgv, 
119              integer *maxdgv, 
120              integer *iordru, 
121              integer *iordrv, 
122              doublereal *xmaxju, 
123              doublereal *xmaxjv, 
124              doublereal *patjac, 
125              doublereal *epmscut, 
126              doublereal *vecerr, 
127              doublereal *erreur, 
128              integer *newdgu, 
129              integer *newdgv);
130
131 static
132 int mma2moy_(integer *ndgumx, 
133              integer *ndgvmx, 
134              integer *ndimen, 
135              integer *mindgu, 
136              integer *maxdgu, 
137              integer *mindgv, 
138              integer *maxdgv, 
139              integer *iordru, 
140              integer *iordrv, 
141              doublereal *patjac, 
142              doublereal *errmoy);
143
144 static
145 int mma2ds2_(integer *ndimen, 
146              doublereal *uintfn, 
147              doublereal *vintfn, 
148              const AdvApp2Var_EvaluatorFunc2Var& foncnp,
149              integer *nbpntu, 
150              integer *nbpntv, 
151              doublereal *urootb, 
152              doublereal *vrootb, 
153              integer *iiuouv, 
154              doublereal *sosotb, 
155              doublereal *disotb, 
156              doublereal *soditb, 
157              doublereal *diditb, 
158              doublereal *fpntab, 
159              doublereal *ttable, 
160              integer *iercod);
161
162
163
164
165 static
166 int mma1fdi_(integer *ndimen, 
167              doublereal *uvfonc, 
168              const AdvApp2Var_EvaluatorFunc2Var& foncnp,
169              integer *isofav, 
170              doublereal *tconst, 
171              integer *nbroot, 
172              doublereal *ttable, 
173              integer *iordre, 
174              integer *ideriv, 
175              doublereal *fpntab, 
176              doublereal *somtab, 
177              doublereal *diftab, 
178              doublereal *contr1,
179              doublereal *contr2, 
180              integer *iercod);
181
182 static
183 int mma1cdi_(integer *ndimen, 
184              integer *nbroot, 
185              doublereal *rootlg, 
186              integer *iordre, 
187              doublereal *contr1, 
188              doublereal *contr2, 
189              doublereal *somtab, 
190              doublereal *diftab, 
191              doublereal *fpntab, 
192              doublereal *hermit, 
193              integer *iercod);
194 static
195 int mma1jak_(integer *ndimen, 
196              integer *nbroot, 
197              integer *iordre,
198              integer *ndgjac, 
199              doublereal *somtab, 
200              doublereal *diftab, 
201              doublereal *cgauss, 
202              doublereal *crvjac, 
203              integer *iercod);
204 static
205 int mma1cnt_(integer *ndimen, 
206              integer *iordre, 
207              doublereal *contr1, 
208              doublereal *contr2, 
209              doublereal *hermit, 
210              integer *ndgjac, 
211              doublereal *crvjac);
212
213 static
214 int mma1fer_(integer *ndimen, 
215              integer *nbsesp, 
216              integer *ndimse, 
217              integer *iordre, 
218              integer *ndgjac, 
219              doublereal *crvjac, 
220              integer *ncflim, 
221              doublereal *epsapr, 
222              doublereal *ycvmax, 
223              doublereal *errmax, 
224              doublereal *errmoy, 
225              integer *ncoeff, 
226              integer *iercod);
227
228 static
229 int mma1noc_(doublereal *dfuvin, 
230              integer *ndimen, 
231              integer *iordre, 
232              doublereal *cntrin, 
233              doublereal *duvout, 
234              integer *isofav, 
235              integer *ideriv, 
236              doublereal *cntout);
237
238
239 static
240   int mmmapcoe_(integer *ndim, 
241                 integer *ndgjac, 
242                 integer *iordre, 
243                 integer *nbpnts, 
244                 doublereal *somtab, 
245                 doublereal *diftab, 
246                 doublereal *gsstab, 
247                 doublereal *crvjac);
248
249 static
250   int mmaperm_(integer *ncofmx, 
251                integer *ndim, 
252                integer *ncoeff, 
253                integer *iordre, 
254                doublereal *crvjac, 
255                integer *ncfnew, 
256                doublereal *errmoy);
257
258
259 #define mmapgss_1 mmapgss_
260 #define mmapgs0_1 mmapgs0_
261 #define mmapgs1_1 mmapgs1_
262 #define mmapgs2_1 mmapgs2_
263
264 //=======================================================================
265 //function : mma1cdi_
266 //purpose  : 
267 //=======================================================================
268 int mma1cdi_(integer *ndimen, 
269              integer *nbroot, 
270              doublereal *rootlg, 
271              integer *iordre, 
272              doublereal *contr1, 
273              doublereal *contr2, 
274              doublereal *somtab, 
275              doublereal *diftab, 
276              doublereal *fpntab, 
277              doublereal *hermit, 
278              integer *iercod)
279 {
280   static integer c__1 = 1;
281
282   /* System generated locals */
283   integer contr1_dim1, contr1_offset, contr2_dim1, contr2_offset, 
284   somtab_dim1, somtab_offset, diftab_dim1, diftab_offset, 
285   fpntab_dim1, fpntab_offset, hermit_dim1, hermit_offset, i__1, 
286   i__2, i__3;
287
288   /* Local variables */
289   static integer nroo2, ncfhe, nd, ii, kk;
290   static integer ibb, kkm, kkp;
291   static doublereal bid1, bid2, bid3;
292
293 /* ********************************************************************** 
294 */
295 /*     FUNCTION : */
296 /*     ---------- */
297 /*     Discretisation on the parameters of interpolation polynomes */
298 /*     constraints of order IORDRE. */
299
300 /*     KEYWORDS : */
301 /*     ----------- */
302 /*     ALL, AB_SPECIFI::CONTRAINTE&, DISCRETISATION, &POINT */
303
304 /*     INPUT ARGUMENTS : */
305 /*     ------------------ */
306 /*     NDIMEN: Space dimension. */
307 /*     NBROOT: Number of INTERNAL discretisation parameters. */
308 /*             It is also the root number Legendre polynome where */
309 /*             the discretization is performed. */
310 /*     ROOTLG: Table of discretization parameters ON (-1,1). */
311 /*     IORDRE: Order of constraint imposed to the extremities of the iso. */
312 /*             = 0, the extremities of the iso are calculated */
313 /*             = 1, additionally, the 1st derivative in the direction */
314 /*                  of the iso is calculated. */
315 /*             = 2, additionally, the 2nd derivative in the direction */
316 /*                  of the iso is calculated. */
317 /*     CONTR1: Contains, if IORDRE>=0, values IORDRE+1 in TTABLE(0) 
318 */
319 /*             (1st extremity) of derivatives of F(Uc,Ve) or F(Ue,Vc), */
320 /*             see below. */
321 /*     CONTR2: Contains, if IORDRE>=0, values IORDRE+1 in */
322 /*             TTABLE(NBROOT+1) (2nd extremity) of: */
323 /*              If ISOFAV=1, derived of order IDERIV by U, derived */
324 /*             ordre 0 to IORDRE by V of F(Uc,Ve) or Uc=TCONST */
325 /*             (fixed iso value) and Ve is the fixed extremity. */
326 /*               If  ISOFAV=2, derivative of order IDERIV by V, derivative */
327 /*             of order 0 to IORDRE by U of F(Ue,Vc) or Vc=TCONST */
328 /*             (fixed iso value) and Ue is the fixed extremity. */
329
330 /*     SOMTAB: Table of NBROOT/2 sums of 2 index points */
331 /*             NBROOT-II+1 and II, for II = 1, NBROOT/2. */
332 /*     DIFTAB: Table of NBROOT/2 differences of 2 index points */
333 /*             NBROOT-II+1 and II, for II = 1, NBROOT/2. */
334
335 /*     OUTPUT ARGUMENTS : */
336 /*     ------------------- */
337 /*     SOMTAB: Table of NBROOT/2 sums of 2 index points */
338 /*             NBROOT-II+1 and II, for II = 1, NBROOT/2 */
339 /*     DIFTAB: Table of  NBROOT/2 differences of 2 index points */
340 /*             NBROOT-II+1 and II, for II = 1, NBROOT/2 */
341 /*     FPNTAB: Auxiliary table. */
342 /*     HERMIT: Table of coeff. 2*(IORDRE+1) Hermite polynoms */
343 /*             of degree 2*IORDRE+1. */
344 /*     IERCOD: Error code, */
345 /*             = 0, Everythig is OK */
346 /*             = 1, The value of IORDRE is out of (0,2) */
347 /*     COMMON USED   : */
348 /*     ---------------- */
349
350 /*     REFERENCES CALLED   : */
351 /*     ----------------------- */
352
353 /*     DESCRIPTION/NOTES/LIMITATIONS : */
354 /*     ----------------------------------- */
355 /*     The results of discretization are arranged in 2 tables */
356 /*     SOMTAB and DIFTAB to earn time during the */
357 /*     calculation of coefficients of the approximation curve. */
358
359 /*     If NBROOT is uneven in SOMTAB(0,*) and DIFTAB(0,*) one stores */
360 /*     the values of the median root of Legendre (0.D0 in (-1,1)). */
361
362 /* ********************************************************************** 
363 */
364
365 /*   Name of the routine */
366
367
368     /* Parameter adjustments */
369     diftab_dim1 = *nbroot / 2 + 1;
370     diftab_offset = diftab_dim1;
371     diftab -= diftab_offset;
372     somtab_dim1 = *nbroot / 2 + 1;
373     somtab_offset = somtab_dim1;
374     somtab -= somtab_offset;
375     --rootlg;
376     hermit_dim1 = (*iordre << 1) + 2;
377     hermit_offset = hermit_dim1;
378     hermit -= hermit_offset;
379     fpntab_dim1 = *nbroot;
380     fpntab_offset = fpntab_dim1 + 1;
381     fpntab -= fpntab_offset;
382     contr2_dim1 = *ndimen;
383     contr2_offset = contr2_dim1 + 1;
384     contr2 -= contr2_offset;
385     contr1_dim1 = *ndimen;
386     contr1_offset = contr1_dim1 + 1;
387     contr1 -= contr1_offset;
388
389     /* Function Body */
390     ibb = AdvApp2Var_SysBase::mnfndeb_();
391     if (ibb >= 3) {
392         AdvApp2Var_SysBase::mgenmsg_("MMA1CDI", 7L);
393     }
394     *iercod = 0;
395
396 /* --- Recuperate 2*(IORDRE+1) coeff of 2*(IORDRE+1) of Hermite polynom --- 
397 */
398
399     AdvApp2Var_ApproxF2var::mma1her_(iordre, &hermit[hermit_offset], iercod);
400     if (*iercod > 0) {
401         goto L9100;
402     }
403
404 /* ------------------- Discretization of Hermite polynoms ----------- 
405 */
406
407     ncfhe = (*iordre + 1) << 1;
408     i__1 = ncfhe;
409     for (ii = 1; ii <= i__1; ++ii) {
410         i__2 = *nbroot;
411         for (kk = 1; kk <= i__2; ++kk) {
412             AdvApp2Var_MathBase::mmmpocur_(&ncfhe, &c__1, &ncfhe, &hermit[ii * hermit_dim1], &
413                     rootlg[kk], &fpntab[kk + ii * fpntab_dim1]);
414 /* L200: */
415         }
416 /* L100: */
417     }
418
419 /* ---- Discretizations of boundary polynoms are taken ---- 
420 */
421
422     nroo2 = *nbroot / 2;
423     i__1 = *ndimen;
424     for (nd = 1; nd <= i__1; ++nd) {
425         i__2 = *iordre + 1;
426         for (ii = 1; ii <= i__2; ++ii) {
427             bid1 = contr1[nd + ii * contr1_dim1];
428             bid2 = contr2[nd + ii * contr2_dim1];
429             i__3 = nroo2;
430             for (kk = 1; kk <= i__3; ++kk) {
431                 kkm = nroo2 - kk + 1;
432                 bid3 = bid1 * fpntab[kkm + ((ii << 1) - 1) * fpntab_dim1] + 
433                         bid2 * fpntab[kkm + (ii << 1) * fpntab_dim1];
434                 somtab[kk + nd * somtab_dim1] -= bid3;
435                 diftab[kk + nd * diftab_dim1] += bid3;
436 /* L500: */
437             }
438             i__3 = nroo2;
439             for (kk = 1; kk <= i__3; ++kk) {
440                 kkp = (*nbroot + 1) / 2 + kk;
441                 bid3 = bid1 * fpntab[kkp + ((ii << 1) - 1) * fpntab_dim1] + 
442                         bid2 * fpntab[kkp + (ii << 1) * fpntab_dim1];
443                 somtab[kk + nd * somtab_dim1] -= bid3;
444                 diftab[kk + nd * diftab_dim1] -= bid3;
445 /* L600: */
446             }
447 /* L400: */
448         }
449 /* L300: */
450     }
451
452 /* ------------ Cas when discretization is done on the roots of a  ----------- 
453 */
454 /* ---------- Legendre polynom of uneven degree, 0 is root -------- 
455 */
456
457     if (*nbroot % 2 == 1) {
458         i__1 = *ndimen;
459         for (nd = 1; nd <= i__1; ++nd) {
460             i__2 = *iordre + 1;
461             for (ii = 1; ii <= i__2; ++ii) {
462                 bid3 = fpntab[nroo2 + 1 + ((ii << 1) - 1) * fpntab_dim1] * 
463                         contr1[nd + ii * contr1_dim1] + fpntab[nroo2 + 1 + (
464                         ii << 1) * fpntab_dim1] * contr2[nd + ii * 
465                         contr2_dim1];
466 /* L800: */
467             }
468             somtab[nd * somtab_dim1] -= bid3;
469             diftab[nd * diftab_dim1] -= bid3;
470 /* L700: */
471         }
472     }
473
474     goto L9999;
475
476 /* ------------------------------ The End ------------------------------- 
477 */
478 /* --> IORDRE is not in the authorized zone. */
479 L9100:
480     *iercod = 1;
481     goto L9999;
482
483 L9999:
484     if (ibb >= 3) {
485         AdvApp2Var_SysBase::mgsomsg_("MMA1CDI", 7L);
486     }
487     return 0;
488 } /* mma1cdi_ */
489
490 //=======================================================================
491 //function : mma1cnt_
492 //purpose  : 
493 //=======================================================================
494 int mma1cnt_(integer *ndimen, 
495              integer *iordre, 
496              doublereal *contr1, 
497              doublereal *contr2, 
498              doublereal *hermit, 
499              integer *ndgjac, 
500              doublereal *crvjac)
501 {
502   /* System generated locals */
503   integer contr1_dim1, contr1_offset, contr2_dim1, contr2_offset, 
504   hermit_dim1, hermit_offset, crvjac_dim1, crvjac_offset, i__1, 
505   i__2, i__3;
506
507   /* Local variables */
508   static integer nd, ii, jj, ibb;
509   static doublereal bid;
510
511
512   /* ***********************************************************************
513    */
514   
515   /*     FUNCTION : */
516   /*     ---------- */
517   /*     Add constraint to polynom. */
518   
519   /*     MOTS CLES : */
520   /*     ----------- */
521   /*     ALL,AB_SPECIFI::COURE&,APPROXIMATION,ADDITION,&CONSTRAINT */
522   
523   /*     INPUT ARGUMENTS : */
524   /*     -------------------- */
525   /*     NDIMEN: Dimension of the space */
526   /*     IORDRE: Order of constraint. */
527   /*     CONTR1: pt of constraint in -1, from order 0 to IORDRE. */
528   /*     CONTR2: Pt of constraint in +1, from order 0 to IORDRE. */
529   /*     HERMIT: Table of Hermit polynoms of order IORDRE. */
530   /*     CRVJAV: Curve of approximation in Jacobi base. */
531   
532   /*     OUTPUT ARGUMENTS : */
533   /*     --------------------- */
534   /*     CRVJAV: Curve of approximation in Jacobi base */
535   /*             to which the polynom of interpolation of constraints is added. */
536   
537   /*     COMMON USED : */
538   /*     ------------------ */
539   
540   
541   /*     REFERENCES CALLED : */
542   /*     --------------------- */
543   
544   
545 /*     DESCRIPTION/NOTES/LIMITATIONS : */
546 /*     ----------------------------------- */
547
548 /* > */
549 /* ***********************************************************************
550  */
551 /*                            DECLARATIONS */
552 /* ***********************************************************************
553  */
554 /*   Name of the routine */
555
556 /* ***********************************************************************
557  */
558 /*                         INITIALISATIONS */
559 /* ***********************************************************************
560  */
561
562     /* Parameter adjustments */
563   hermit_dim1 = (*iordre << 1) + 2;
564   hermit_offset = hermit_dim1;
565   hermit -= hermit_offset;
566   contr2_dim1 = *ndimen;
567   contr2_offset = contr2_dim1 + 1;
568   contr2 -= contr2_offset;
569   contr1_dim1 = *ndimen;
570   contr1_offset = contr1_dim1 + 1;
571   contr1 -= contr1_offset;
572   crvjac_dim1 = *ndgjac + 1;
573   crvjac_offset = crvjac_dim1;
574   crvjac -= crvjac_offset;
575   
576   /* Function Body */
577   ibb = AdvApp2Var_SysBase::mnfndeb_();
578   if (ibb >= 3) {
579     AdvApp2Var_SysBase::mgenmsg_("MMA1CNT", 7L);
580     }
581
582 /* ***********************************************************************
583  */
584 /*                            Processing */
585 /* ***********************************************************************
586  */
587
588   i__1 = *ndimen;
589   for (nd = 1; nd <= i__1; ++nd) {
590     i__2 = (*iordre << 1) + 1;
591     for (ii = 0; ii <= i__2; ++ii) {
592       bid = 0.;
593       i__3 = *iordre + 1;
594       for (jj = 1; jj <= i__3; ++jj) {
595         bid = bid + contr1[nd + jj * contr1_dim1] *
596           hermit[ii + ((jj  << 1) - 1) * hermit_dim1] + 
597             contr2[nd + jj * contr2_dim1] * hermit[ii + (jj << 1) * hermit_dim1];
598         /* L300: */
599       }
600       crvjac[ii + nd * crvjac_dim1] = bid;
601       /* L200: */
602     }
603     /* L100: */
604   }
605
606 /* ***********************************************************************
607  */
608 /*                   RETURN CALLING PROGRAM */
609 /* ***********************************************************************
610  */
611
612     if (ibb >= 3) {
613         AdvApp2Var_SysBase::mgsomsg_("MMA1CNT", 7L);
614     }
615
616   return 0 ;
617 } /* mma1cnt_ */
618
619 //=======================================================================
620 //function : mma1fdi_
621 //purpose  : 
622 //=======================================================================
623 int mma1fdi_(integer *ndimen, 
624              doublereal *uvfonc, 
625              const AdvApp2Var_EvaluatorFunc2Var& foncnp,
626              integer *isofav, 
627              doublereal *tconst, 
628              integer *nbroot, 
629              doublereal *ttable, 
630              integer *iordre, 
631              integer *ideriv, 
632              doublereal *fpntab, 
633              doublereal *somtab, 
634              doublereal *diftab, 
635              doublereal *contr1,
636              doublereal *contr2, 
637              integer *iercod)
638 {
639   /* System generated locals */
640   integer fpntab_dim1, somtab_dim1, somtab_offset, diftab_dim1, 
641   diftab_offset, contr1_dim1, contr1_offset, contr2_dim1, 
642   contr2_offset, i__1, i__2;
643   doublereal d__1;
644
645   /* Local variables */
646   static integer ideb, ifin, nroo2, ideru, iderv;
647   static doublereal renor;
648   static integer ii, nd, ibb, iim, nbp, iip;
649   static doublereal bid1, bid2;
650
651 /* ********************************************************************** 
652 */
653
654 /*     FUNCTION : */
655 /*     ---------- */
656 /*     DiscretiZation of a non-polynomial function F(U,V) or of */
657 /*     its derivative with fixed isoparameter. */
658
659 /*     KEYWORDS : */
660 /*     ----------- */
661 /*     ALL, AB_SPECIFI::FONCTION&, DISCRETISATION, &POINT */
662
663 /*     INPUT ARGUMENTS : */
664 /*     ------------------ */
665 /*     NDIMEN: Space dimension. */
666 /*     UVFONC: Limits of the path of definition by U and by V of the approximated function */
667 /*     FONCNP: The NAME of the non-polynomial function to be approximated */
668 /*             (external program). */
669 /*     ISOFAV: Fixed isoparameter for the discretization; */
670 /*             = 1, discretization with fixed U and variable V. */
671 /*             = 2, discretization with fixed V and variable U. */
672 /*     TCONST: Iso value is also fixed. */
673 /*     NBROOT: Number of INTERNAL discretization parameters. */
674 /*             (if there are constraints, 2 extremities should be added). 
675 */
676 /*             This is also the root number of the Legendre polynom where */
677 /*             the discretization is done. */
678 /*     TTABLE: Table of discretization parameters and of 2 extremities */
679 /*             (Respectively (-1, NBROOT Legendre roots,1) */
680 /*             reframed within the adequate interval. */
681 /*     IORDRE: Order of constraint imposed on the extremities of the iso. */
682 /*             (If Iso-U, it is necessary to calculate the derivatives by V and vice */
683 /*             versa). */
684 /*             = 0, the extremities of the iso are calculated. */
685 /*             = 1, additionally the 1st derivative in the direction of the iso is calculated */
686 /*             = 2, additionally the 2nd derivative in the direction of the iso is calculated */
687 /*     IDERIV: Order of derivative transversal to fixed iso (If Iso-U=Uc */
688 /*             is fixed, the derivative of order IDERIV is discretized by U of */
689 /*             F(Uc,v). Same if iso-V is fixed). */
690 /*             Varies from 0 (positioning) to 2 (2nd derivative). */
691
692 /*     OUTPUT ARGUMENTS : */
693 /*     ------------------- */
694 /*     FPNTAB: Auxiliary table. 
695        SOMTAB: Table of NBROOT/2 sums of 2 index points */
696 /*             NBROOT-II+1 and II, for II = 1, NBROOT/2 */
697 /*     DIFTAB: Table of  NBROOT/2 differences of 2 index points */
698 /*             NBROOT-II+1 and II, for II = 1, NBROOT/2 */
699 /*     CONTR1: Contains, if IORDRE>=0, values IORDRE+1 in TTABLE(0) 
700 */
701 /*             (1st extremity) of derivatives of F(Uc,Ve) or F(Ue,Vc), */
702 /*             see below. */
703 /*     CONTR2: Contains, if IORDRE>=0, values IORDRE+1 in */
704 /*             TTABLE(NBROOT+1) (2nd extremity) of: */
705 /*              If ISOFAV=1, derived of order IDERIV by U, derived */
706 /*             ordre 0 to IORDRE by V of F(Uc,Ve) or Uc=TCONST */
707 /*             (fixed iso value) and Ve is the fixed extremity. */
708 /*               If  ISOFAV=2, derivative of order IDERIV by V, derivative */
709 /*             of order 0 to IORDRE by U of F(Ue,Vc) or Vc=TCONST */
710 /*             (fixed iso value) and Ue is the fixed extremity. */
711 /*     IERCOD: Error code > 100; Pb in  evaluation of FONCNP, */
712 /*             the returned error code is equal to error code of FONCNP + 100. */
713
714 /*     COMMONS USED   : */
715 /*     ---------------- */
716
717 /*     REFERENCES CALLED   : */
718 /*     ----------------------- */
719
720 /*     DESCRIPTION/NOTES/LIMITATIONS : */
721 /*     ----------------------------------- */
722 /*     The results of discretization are arranged in 2 tables */
723 /*     SOMTAB and DIFTAB to earn time during the */
724 /*     calculation of coefficients of the approximation curve. */
725
726 /*     If NBROOT is uneven in SOMTAB(0,*) and DIFTAB(0,*) one stores */
727 /*     the values of the median root of Legendre (0.D0 in (-1,1)). */
728
729 /*     Function F(u,v) defined in UVFONC is reparameterized in */
730 /*     (-1,1)x(-1,1). Then 1st and 2nd derivatives are renormalized. */
731
732 /* > */
733 /* ********************************************************************** 
734 */
735
736 /*   Name of the routine */
737
738
739     /* Parameter adjustments */
740     uvfonc -= 3;
741     diftab_dim1 = *nbroot / 2 + 1;
742     diftab_offset = diftab_dim1;
743     diftab -= diftab_offset;
744     somtab_dim1 = *nbroot / 2 + 1;
745     somtab_offset = somtab_dim1;
746     somtab -= somtab_offset;
747     fpntab_dim1 = *ndimen;
748     --fpntab;
749     contr2_dim1 = *ndimen;
750     contr2_offset = contr2_dim1 + 1;
751     contr2 -= contr2_offset;
752     contr1_dim1 = *ndimen;
753     contr1_offset = contr1_dim1 + 1;
754     contr1 -= contr1_offset;
755
756     /* Function Body */
757     ibb = AdvApp2Var_SysBase::mnfndeb_();
758     if (ibb >= 3) {
759         AdvApp2Var_SysBase::mgenmsg_("MMA1FDI", 7L);
760     }
761     *iercod = 0;
762
763 /* --------------- Definition of the nb of points to calculate -------------- 
764 */
765 /* --> If constraints, the limits are also taken */
766     if (*iordre >= 0) {
767         ideb = 0;
768         ifin = *nbroot + 1;
769 /* --> Otherwise, only Legendre roots (reframed) are used
770 . */
771     } else {
772         ideb = 1;
773         ifin = *nbroot;
774     }
775 /* --> Nb of point to calculate. */
776     nbp = ifin - ideb + 1;
777     nroo2 = *nbroot / 2;
778
779 /* --------------- Determination of the order of global derivation -------- 
780 */
781 /* --> ISOFAV takes only values 1 or 2. */
782 /*    if Iso-U, derive by U of order IDERIV */
783     if (*isofav == 1) {
784         ideru = *ideriv;
785         iderv = 0;
786         d__1 = (uvfonc[4] - uvfonc[3]) / 2.;
787         renor = AdvApp2Var_MathBase::pow__di(&d__1, ideriv);
788 /*    if Iso-V, derive by V of order IDERIV */
789     } else {
790         ideru = 0;
791         iderv = *ideriv;
792         d__1 = (uvfonc[6] - uvfonc[5]) / 2.;
793         renor = AdvApp2Var_MathBase::pow__di(&d__1, ideriv);
794     }
795
796 /* ----------- Discretization on roots of the  --------------- 
797 */
798 /* ---------------------- Legendre polynom of degree NBROOT ------------------- 
799 */
800
801     (*const_cast <AdvApp2Var_EvaluatorFunc2Var*> (&foncnp)).Evaluate (ndimen, 
802               &uvfonc[3], 
803               &uvfonc[5], 
804               isofav, 
805               tconst, 
806               &nbp, 
807               &ttable[ideb], 
808               &ideru, 
809               &iderv, 
810               &fpntab[ideb * fpntab_dim1 + 1], 
811               iercod);
812     if (*iercod > 0) {
813         goto L9999;
814     }
815     i__1 = *ndimen;
816     for (nd = 1; nd <= i__1; ++nd) {
817         i__2 = nroo2;
818         for (ii = 1; ii <= i__2; ++ii) {
819             iip = (*nbroot + 1) / 2 + ii;
820             iim = nroo2 - ii + 1;
821             bid1 = fpntab[nd + iim * fpntab_dim1];
822             bid2 = fpntab[nd + iip * fpntab_dim1];
823             somtab[ii + nd * somtab_dim1] = renor * (bid2 + bid1);
824             diftab[ii + nd * diftab_dim1] = renor * (bid2 - bid1);
825 /* L200: */
826         }
827 /* L100: */
828     }
829
830 /* ------------ Case when discretisation is done on roots of a ---- 
831 */
832 /* ---------- Legendre polynom of uneven degree, 0 is root -------- 
833 */
834
835     if (*nbroot % 2 == 1) {
836         i__1 = *ndimen;
837         for (nd = 1; nd <= i__1; ++nd) {
838             somtab[nd * somtab_dim1] = renor * fpntab[nd + (nroo2 + 1) * 
839                     fpntab_dim1];
840             diftab[nd * diftab_dim1] = renor * fpntab[nd + (nroo2 + 1) * 
841                     fpntab_dim1];
842 /* L300: */
843         }
844     } else {
845         i__1 = *ndimen;
846         for (nd = 1; nd <= i__1; ++nd) {
847             somtab[nd * somtab_dim1] = 0.;
848             diftab[nd * diftab_dim1] = 0.;
849         }
850     }
851
852
853 /* --------------------- Take into account constraints ---------------- 
854 */
855
856     if (*iordre >= 0) {
857 /* --> Recover already calculated extremities. */
858         i__1 = *ndimen;
859         for (nd = 1; nd <= i__1; ++nd) {
860             contr1[nd + contr1_dim1] = renor * fpntab[nd];
861             contr2[nd + contr2_dim1] = renor * fpntab[nd + (*nbroot + 1) * 
862                     fpntab_dim1];
863 /* L400: */
864         }
865 /* --> Nb of points to calculate/call to FONCNP */
866         nbp = 1;
867 /*    If Iso-U, derive by V till order IORDRE */
868         if (*isofav == 1) {
869 /* --> Factor of normalisation 1st derivative. */
870             bid1 = (uvfonc[6] - uvfonc[5]) / 2.;
871             i__1 = *iordre;
872             for (iderv = 1; iderv <= i__1; ++iderv) {
873                 (*const_cast <AdvApp2Var_EvaluatorFunc2Var*> (&foncnp)).Evaluate (
874             ndimen, &uvfonc[3], &uvfonc[5], isofav, tconst, &
875                         nbp, ttable, &ideru, &iderv, &contr1[(iderv + 1) * 
876                         contr1_dim1 + 1], iercod);
877                 if (*iercod > 0) {
878                     goto L9999;
879                 }
880 /* L500: */
881             }
882             i__1 = *iordre;
883             for (iderv = 1; iderv <= i__1; ++iderv) {
884                 (*const_cast <AdvApp2Var_EvaluatorFunc2Var*> (&foncnp)).Evaluate (
885             ndimen, &uvfonc[3], &uvfonc[5], isofav, tconst, &
886                         nbp, &ttable[*nbroot + 1], &ideru, &iderv, &contr2[(
887                         iderv + 1) * contr2_dim1 + 1], iercod);
888                 if (*iercod > 0) {
889                     goto L9999;
890                 }
891 /* L510: */
892             }
893 /*    If Iso-V, derive by U till order IORDRE */
894         } else {
895 /* --> Factor of normalization  1st derivative. */
896             bid1 = (uvfonc[4] - uvfonc[3]) / 2.;
897             i__1 = *iordre;
898             for (ideru = 1; ideru <= i__1; ++ideru) {
899                 (*const_cast <AdvApp2Var_EvaluatorFunc2Var*> (&foncnp)).Evaluate (
900             ndimen, &uvfonc[3], &uvfonc[5], isofav, tconst, &
901                         nbp, ttable, &ideru, &iderv, &contr1[(ideru + 1) * 
902                         contr1_dim1 + 1], iercod);
903                 if (*iercod > 0) {
904                     goto L9999;
905                 }
906 /* L600: */
907             }
908             i__1 = *iordre;
909             for (ideru = 1; ideru <= i__1; ++ideru) {
910                 (*const_cast <AdvApp2Var_EvaluatorFunc2Var*> (&foncnp)).Evaluate (
911             ndimen, &uvfonc[3], &uvfonc[5], isofav, tconst, &
912                         nbp, &ttable[*nbroot + 1], &ideru, &iderv, &contr2[(
913                         ideru + 1) * contr2_dim1 + 1], iercod);
914                 if (*iercod > 0) {
915                     goto L9999;
916                 }
917 /* L610: */
918             }
919         }
920
921 /* ------------------------- Normalization of derivatives -------------
922 ---- */
923 /* (The function is redefined on (-1,1)*(-1,1)) */
924         bid2 = renor;
925         i__1 = *iordre;
926         for (ii = 1; ii <= i__1; ++ii) {
927             bid2 = bid1 * bid2;
928             i__2 = *ndimen;
929             for (nd = 1; nd <= i__2; ++nd) {
930                 contr1[nd + (ii + 1) * contr1_dim1] *= bid2;
931                 contr2[nd + (ii + 1) * contr2_dim1] *= bid2;
932 /* L710: */
933             }
934 /* L700: */
935         }
936     }
937
938 /* ------------------------------ The end ------------------------------- 
939 */
940
941 L9999:
942     if (*iercod > 0) {
943         *iercod += 100;
944         AdvApp2Var_SysBase::maermsg_("MMA1FDI", iercod, 7L);
945     }
946     if (ibb >= 3) {
947         AdvApp2Var_SysBase::mgsomsg_("MMA1FDI", 7L);
948     }
949     return 0;
950 } /* mma1fdi_ */
951
952 //=======================================================================
953 //function : mma1fer_
954 //purpose  : 
955 //=======================================================================
956 int mma1fer_(integer *,//ndimen, 
957              integer *nbsesp, 
958              integer *ndimse, 
959              integer *iordre, 
960              integer *ndgjac, 
961              doublereal *crvjac, 
962              integer *ncflim, 
963              doublereal *epsapr, 
964              doublereal *ycvmax, 
965              doublereal *errmax, 
966              doublereal *errmoy, 
967              integer *ncoeff, 
968              integer *iercod)
969 {
970   /* System generated locals */
971   integer crvjac_dim1, crvjac_offset, i__1, i__2;
972
973   /* Local variables */
974   static integer idim, ncfja, ncfnw, ndses, ii, kk, ibb, ier;
975   static integer nbr0;
976
977
978 /* ***********************************************************************
979  */
980
981 /*     FUNCTION : */
982 /*     ---------- */
983 /*     Calculate the degree and the errors of approximation of a border. */
984
985 /*     KEYWORDS : */
986 /*     ----------- */
987 /*      TOUS,AB_SPECIFI :: COURBE&,TRONCATURE, &PRECISION */
988
989 /*     INPUT ARGUMENTS : */
990 /*     -------------------- */
991
992 /*     NDIMEN: Total Dimension of the space (sum of dimensions of sub-spaces) */
993 /*     NBSESP: Number of "independent" sub-spaces. */
994 /*     NDIMSE: Table of dimensions of sub-spaces. */
995 /*     IORDRE: Order of constraint at the extremities of the border */
996 /*              -1 = no constraints, */
997 /*               0 = constraints of passage to limits (i.e. C0), */
998 /*               1 = C0 + constraintes of 1st derivatives (i.e. C1), */
999 /*               2 = C1 + constraintes of 2nd derivatives (i.e. C2). */
1000 /*     NDGJAC: Degree of development in series to use for the calculation  
1001 /*             in the base of Jacobi. */
1002 /*     CRVJAC: Table of coeff. of the curve of approximation in the */
1003 /*             base of Jacobi. */
1004 /*     NCFLIM: Max number of coeff of the polynomial curve */
1005 /*             of approximation (should be above or equal to */
1006 /*             2*IORDRE+2 and below or equal to 50). */
1007 /*     EPSAPR: Table of errors of approximations that cannot be passed, */
1008 /*             sub-space by sub-space. */
1009
1010 /*     OUTPUT ARGUMENTS : */
1011 /*     --------------------- */
1012 /*     YCVMAX: Auxiliary Table. */
1013 /*     ERRMAX: Table of errors (sub-space by sub-space) */
1014 /*             MAXIMUM made in the approximation of FONCNP by */
1015 /*             COURBE. */
1016 /*     ERRMOY: Table of errors (sub-space by sub-space) */
1017 /*             AVERAGE made in the approximation of FONCNP by */
1018 /*             COURBE. */
1019 /*     NCOEFF: Number of significative coeffs. of the calculated "curve". */
1020 /*     IERCOD: Error code */
1021 /*             = 0, ok, */
1022 /*             =-1, warning, required tolerance can't be */
1023 /*                  met with coefficients NFCLIM. */
1024 /*             = 1, order of constraints (IORDRE) is not within authorised values */
1025 /*                  
1026
1027 /*     COMMONS USED : */
1028 /*     ------------------ */
1029
1030 /*     REFERENCES CALLED : */
1031 /*     --------------------- */
1032
1033 /*     DESCRIPTION/NOTES/LIMITATIONS : */
1034 /*     ----------------------------------- */
1035 /* > */
1036 /* ********************************************************************** 
1037 */
1038
1039 /*  Name of the routine */
1040
1041
1042     /* Parameter adjustments */
1043     --ycvmax;
1044     --errmoy;
1045     --errmax;
1046     --epsapr;
1047     --ndimse;
1048     crvjac_dim1 = *ndgjac + 1;
1049     crvjac_offset = crvjac_dim1;
1050     crvjac -= crvjac_offset;
1051
1052     /* Function Body */
1053     ibb = AdvApp2Var_SysBase::mnfndeb_();
1054     if (ibb >= 3) {
1055         AdvApp2Var_SysBase::mgenmsg_("MMA1FER", 7L);
1056     }
1057     *iercod = 0;
1058     idim = 1;
1059     *ncoeff = 0;
1060     ncfja = *ndgjac + 1;
1061
1062 /* ------------ Calculate the degree of the curve and of the Max error -------- 
1063 */
1064 /* -------------- of approximation for all sub-spaces -------- 
1065 */
1066
1067     i__1 = *nbsesp;
1068     for (ii = 1; ii <= i__1; ++ii) {
1069         ndses = ndimse[ii];
1070
1071 /* ------------ cutting of coeff. and calculation of Max error -------
1072 ---- */
1073
1074         AdvApp2Var_MathBase::mmtrpjj_(&ncfja, &ndses, &ncfja, &epsapr[ii], iordre, &crvjac[idim * 
1075                 crvjac_dim1], &ycvmax[1], &errmax[ii], &ncfnw);
1076
1077 /* ******************************************************************
1078 **** */
1079 /* ------------- If precision OK, calculate the average error -------
1080 ---- */
1081 /* ******************************************************************
1082 **** */
1083
1084         if (ncfnw <= *ncflim) {
1085             mmaperm_(&ncfja, &ndses, &ncfja, iordre, &crvjac[idim * 
1086                     crvjac_dim1], &ncfnw, &errmoy[ii]);
1087             *ncoeff = advapp_max(ncfnw,*ncoeff);
1088
1089 /* ------------- Set the declined coefficients to 0.D0 -----------
1090 -------- */
1091
1092             nbr0 = *ncflim - ncfnw;
1093             if (nbr0 > 0) {
1094                 i__2 = ndses;
1095                 for (kk = 1; kk <= i__2; ++kk) {
1096                   AdvApp2Var_SysBase::mvriraz_(&nbr0, 
1097                              &crvjac[ncfnw + (idim + kk - 1) * crvjac_dim1]);
1098 /* L200: */
1099                 }
1100             }
1101         } else {
1102
1103 /* **************************************************************
1104 ******** */
1105 /* ------------------- If required precision can't be reached----
1106 -------- */
1107 /* **************************************************************
1108 ******** */
1109
1110             *iercod = -1;
1111
1112 /* ------------------------- calculate the Max error ------------
1113 -------- */
1114
1115             AdvApp2Var_MathBase::mmaperx_(&ncfja, &ndses, &ncfja, iordre, &crvjac[idim * 
1116                     crvjac_dim1], ncflim, &ycvmax[1], &errmax[ii], &ier);
1117             if (ier > 0) {
1118                 goto L9100;
1119             }
1120
1121 /* -------------------- nb of coeff to be returned -------------
1122 -------- */
1123
1124             *ncoeff = *ncflim;
1125
1126 /* ------------------- and calculation of the average error ----
1127 -------- */
1128
1129             mmaperm_(&ncfja, &ndses, &ncfja, iordre, &crvjac[idim * 
1130                     crvjac_dim1], ncflim, &errmoy[ii]);
1131         }
1132         idim += ndses;
1133 /* L100: */
1134     }
1135
1136     goto L9999;
1137
1138 /* ------------------------------ The end ------------------------------- 
1139 */
1140 /* --> The order of constraints is not within autorized values. */
1141 L9100:
1142     *iercod = 1;
1143     goto L9999;
1144
1145 L9999:
1146     if (*iercod != 0) {
1147         AdvApp2Var_SysBase::maermsg_("MMA1FER", iercod, 7L);
1148     }
1149     if (ibb >= 3) {
1150         AdvApp2Var_SysBase::mgsomsg_("MMA1FER", 7L);
1151     }
1152     return 0;
1153 } /* mma1fer_ */
1154
1155
1156 //=======================================================================
1157 //function : mma1her_
1158 //purpose  : 
1159 //=======================================================================
1160 int AdvApp2Var_ApproxF2var::mma1her_(const integer *iordre, 
1161                                      doublereal *hermit, 
1162                                      integer *iercod)
1163 {
1164   /* System generated locals */
1165   integer hermit_dim1, hermit_offset;
1166
1167   /* Local variables */
1168   static integer ibb;
1169
1170
1171
1172 /* ********************************************************************** 
1173 */
1174
1175 /*     FUNCTION : */
1176 /*     ---------- */
1177 /*     Calculate 2*(IORDRE+1) Hermit polynoms of  degree 2*IORDRE+1 */
1178 /*     on (-1,1) */
1179
1180 /*     KEYWORDS : */
1181 /*     ----------- */
1182 /*     ALL, AB_SPECIFI::CONTRAINTE&, INTERPOLATION, &POLYNOME */
1183
1184 /*     INPUT ARGUMENTS : */
1185 /*     ------------------ */
1186 /*     IORDRE: Order of constraint. */
1187 /*      = 0, Polynom of interpolation of order C0 on (-1,1). */
1188 /*      = 1, Polynom of interpolation of order C0 and C1 on (-1,1). */
1189 /*      = 2, Polynom of interpolation of order C0, C1 and C2 on (-1,1). 
1190 */
1191
1192 /*     OUTPUT ARGUMENTS : */
1193 /*     ------------------- */
1194 /*     HERMIT: Table of 2*IORDRE+2 coeff. of each of  2*(IORDRE+1) */
1195 /*             HERMIT polynom. */
1196 /*     IERCOD: Error code, */
1197 /*      = 0, Ok */
1198 /*      = 1, required order of constraint is not managed here. */
1199 /*     COMMONS USED   : */
1200 /*     ---------------- */
1201
1202 /*     REFERENCES CALLED   : */
1203 /*     ----------------------- */
1204
1205 /*     DESCRIPTION/NOTES/LIMITATIONS : */
1206 /*     ----------------------------------- */
1207 /*     The part of HERMIT(*,2*i+j) table where  j=1 or 2 and i=0 to IORDRE, 
1208 /*     contains the coefficients of the polynom of degree 2*IORDRE+1 */
1209 /*     such as ALL values in -1 and in +1 of this polynom and its */
1210 /*     derivatives till order of derivation IORDRE are NULL, */
1211 /*     EXCEPT for the derivative of order i: */
1212 /*        - valued 1 in -1 if j=1 */
1213 /*        - valued 1 in +1 if j=2. */
1214 /* > */
1215 /* ********************************************************************** 
1216 */
1217
1218 /*  Name of the routine */
1219
1220
1221     /* Parameter adjustments */
1222     hermit_dim1 = (*iordre + 1) << 1;
1223     hermit_offset = hermit_dim1 + 1;
1224     hermit -= hermit_offset;
1225
1226     /* Function Body */
1227     ibb = AdvApp2Var_SysBase::mnfndeb_();
1228     if (ibb >= 3) {
1229         AdvApp2Var_SysBase::mgenmsg_("MMA1HER", 7L);
1230     }
1231     *iercod = 0;
1232
1233 /* --- Recover (IORDRE+2) coeff of 2*(IORDRE+1) Hermit polynoms -- 
1234 */
1235
1236     if (*iordre == 0) {
1237         hermit[hermit_dim1 + 1] = .5;
1238         hermit[hermit_dim1 + 2] = -.5;
1239
1240         hermit[(hermit_dim1 << 1) + 1] = .5;
1241         hermit[(hermit_dim1 << 1) + 2] = .5;
1242     } else if (*iordre == 1) {
1243         hermit[hermit_dim1 + 1] = .5;
1244         hermit[hermit_dim1 + 2] = -.75;
1245         hermit[hermit_dim1 + 3] = 0.;
1246         hermit[hermit_dim1 + 4] = .25;
1247
1248         hermit[(hermit_dim1 << 1) + 1] = .5;
1249         hermit[(hermit_dim1 << 1) + 2] = .75;
1250         hermit[(hermit_dim1 << 1) + 3] = 0.;
1251         hermit[(hermit_dim1 << 1) + 4] = -.25;
1252
1253         hermit[hermit_dim1 * 3 + 1] = .25;
1254         hermit[hermit_dim1 * 3 + 2] = -.25;
1255         hermit[hermit_dim1 * 3 + 3] = -.25;
1256         hermit[hermit_dim1 * 3 + 4] = .25;
1257
1258         hermit[(hermit_dim1 << 2) + 1] = -.25;
1259         hermit[(hermit_dim1 << 2) + 2] = -.25;
1260         hermit[(hermit_dim1 << 2) + 3] = .25;
1261         hermit[(hermit_dim1 << 2) + 4] = .25;
1262     } else if (*iordre == 2) {
1263         hermit[hermit_dim1 + 1] = .5;
1264         hermit[hermit_dim1 + 2] = -.9375;
1265         hermit[hermit_dim1 + 3] = 0.;
1266         hermit[hermit_dim1 + 4] = .625;
1267         hermit[hermit_dim1 + 5] = 0.;
1268         hermit[hermit_dim1 + 6] = -.1875;
1269
1270         hermit[(hermit_dim1 << 1) + 1] = .5;
1271         hermit[(hermit_dim1 << 1) + 2] = .9375;
1272         hermit[(hermit_dim1 << 1) + 3] = 0.;
1273         hermit[(hermit_dim1 << 1) + 4] = -.625;
1274         hermit[(hermit_dim1 << 1) + 5] = 0.;
1275         hermit[(hermit_dim1 << 1) + 6] = .1875;
1276
1277         hermit[hermit_dim1 * 3 + 1] = .3125;
1278         hermit[hermit_dim1 * 3 + 2] = -.4375;
1279         hermit[hermit_dim1 * 3 + 3] = -.375;
1280         hermit[hermit_dim1 * 3 + 4] = .625;
1281         hermit[hermit_dim1 * 3 + 5] = .0625;
1282         hermit[hermit_dim1 * 3 + 6] = -.1875;
1283
1284         hermit[(hermit_dim1 << 2) + 1] = -.3125;
1285         hermit[(hermit_dim1 << 2) + 2] = -.4375;
1286         hermit[(hermit_dim1 << 2) + 3] = .375;
1287         hermit[(hermit_dim1 << 2) + 4] = .625;
1288         hermit[(hermit_dim1 << 2) + 5] = -.0625;
1289         hermit[(hermit_dim1 << 2) + 6] = -.1875;
1290
1291         hermit[hermit_dim1 * 5 + 1] = .0625;
1292         hermit[hermit_dim1 * 5 + 2] = -.0625;
1293         hermit[hermit_dim1 * 5 + 3] = -.125;
1294         hermit[hermit_dim1 * 5 + 4] = .125;
1295         hermit[hermit_dim1 * 5 + 5] = .0625;
1296         hermit[hermit_dim1 * 5 + 6] = -.0625;
1297
1298         hermit[hermit_dim1 * 6 + 1] = .0625;
1299         hermit[hermit_dim1 * 6 + 2] = .0625;
1300         hermit[hermit_dim1 * 6 + 3] = -.125;
1301         hermit[hermit_dim1 * 6 + 4] = -.125;
1302         hermit[hermit_dim1 * 6 + 5] = .0625;
1303         hermit[hermit_dim1 * 6 + 6] = .0625;
1304     } else {
1305         *iercod = 1;
1306     }
1307
1308 /* ------------------------------ The End ------------------------------- 
1309 */
1310
1311     AdvApp2Var_SysBase::maermsg_("MMA1HER", iercod, 7L);
1312     if (ibb >= 3) {
1313         AdvApp2Var_SysBase::mgsomsg_("MMA1HER", 7L);
1314     }
1315     return 0;
1316 } /* mma1her_ */
1317 //=======================================================================
1318 //function : mma1jak_
1319 //purpose  : 
1320 //=======================================================================
1321 int mma1jak_(integer *ndimen, 
1322              integer *nbroot, 
1323              integer *iordre,
1324              integer *ndgjac, 
1325              doublereal *somtab, 
1326              doublereal *diftab, 
1327              doublereal *cgauss, 
1328              doublereal *crvjac, 
1329              integer *iercod)
1330 {
1331   /* System generated locals */
1332   integer somtab_dim1, somtab_offset, diftab_dim1, diftab_offset, 
1333   crvjac_dim1, crvjac_offset, cgauss_dim1;
1334
1335   /* Local variables */
1336   static integer ibb;
1337
1338 /* ********************************************************************** 
1339 */
1340
1341 /*     FUNCTION : */
1342 /*     ---------- */
1343 /*     Calculate the curve of approximation of a non-polynomial function */
1344 /*     in the base of Jacobi. */
1345
1346 /*     KEYWORDS : */
1347 /*     ----------- */
1348 /*     FUNCTION,DISCRETISATION,APPROXIMATION,CONSTRAINT,CURVE,JACOBI */
1349
1350 /*     INPUT ARGUMENTS : */
1351 /*     ------------------ */
1352 /*     NDIMEN: Total dimension of the space (sum of dimensions */
1353 /*             of sub-spaces) */
1354 /*     NBROOT: Nb of points of discretization of the iso, extremities not 
1355 /*             included. */
1356 /*     IORDRE: Order of constraint at the extremities of the boundary */
1357 /*              -1 = no constraints, */
1358 /*               0 = constraints of passage of limits (i.e. C0), */
1359 /*               1 = C0 + constraints of 1st derivatives (i.e. C1), */
1360 /*               2 = C1 + constraints of 2nd derivatives (i.e. C2). */
1361 /*     NDGJAC: Degree of development in series to be used for calculation in the  
1362 /*             base of Jacobi. */
1363
1364 /*     OUTPUT ARGUMENTS : */
1365 /*     ------------------- */
1366 /*     CRVJAC : Curve of approximation of FONCNP with (eventually) */
1367 /*              taking into account of constraints at the extremities. */
1368 /*              This curve is of degree NDGJAC. */
1369 /*     IERCOD : Error code : */
1370 /*               0 = All is ok. */
1371 /*              33 = Pb to return data of du block data */
1372 /*                   of coeff. of integration by GAUSS method. */
1373 /*                   by program MMAPPTT. */
1374
1375 /*     COMMONS USED   : */
1376 /*     ---------------- */
1377
1378 /*     REFERENCES CALLED   : */
1379 /*     ----------------------- */
1380 /* > */
1381 /* ********************************************************************** 
1382 */
1383
1384 /*   Name of the routine */
1385
1386     /* Parameter adjustments */
1387     diftab_dim1 = *nbroot / 2 + 1;
1388     diftab_offset = diftab_dim1;
1389     diftab -= diftab_offset;
1390     somtab_dim1 = *nbroot / 2 + 1;
1391     somtab_offset = somtab_dim1;
1392     somtab -= somtab_offset;
1393     crvjac_dim1 = *ndgjac + 1;
1394     crvjac_offset = crvjac_dim1;
1395     crvjac -= crvjac_offset;
1396     cgauss_dim1 = *nbroot / 2 + 1;
1397
1398     /* Function Body */
1399     ibb = AdvApp2Var_SysBase::mnfndeb_();
1400     if (ibb >= 2) {
1401         AdvApp2Var_SysBase::mgenmsg_("MMA1JAK", 7L);
1402     }
1403     *iercod = 0;
1404
1405 /* ----------------- Recover coeffs of integration by Gauss ----------- 
1406 */
1407
1408     AdvApp2Var_ApproxF2var::mmapptt_(ndgjac, nbroot, iordre, cgauss, iercod);
1409     if (*iercod > 0) {
1410         *iercod = 33;
1411         goto L9999;
1412     }
1413
1414 /* --------------- Calculate the curve in the base of Jacobi ----------- 
1415 */
1416
1417     mmmapcoe_(ndimen, ndgjac, iordre, nbroot, &somtab[somtab_offset], &diftab[
1418             diftab_offset], cgauss, &crvjac[crvjac_offset]);
1419
1420 /* ------------------------------ The End ------------------------------- 
1421 */
1422
1423 L9999:
1424     if (*iercod != 0) {
1425         AdvApp2Var_SysBase::maermsg_("MMA1JAK", iercod, 7L);
1426     }
1427     if (ibb >= 2) {
1428         AdvApp2Var_SysBase::mgsomsg_("MMA1JAK", 7L);
1429     }
1430     return 0;
1431 } /* mma1jak_ */
1432
1433 //=======================================================================
1434 //function : mma1noc_
1435 //purpose  : 
1436 //=======================================================================
1437 int mma1noc_(doublereal *dfuvin, 
1438              integer *ndimen, 
1439              integer *iordre, 
1440              doublereal *cntrin, 
1441              doublereal *duvout, 
1442              integer *isofav, 
1443              integer *ideriv, 
1444              doublereal *cntout)
1445 {
1446   /* System generated locals */
1447   integer i__1;
1448   doublereal d__1;
1449
1450   /* Local variables */
1451   static doublereal rider, riord;
1452   static integer nd, ibb;
1453   static doublereal bid;
1454 /* ********************************************************************** 
1455 */
1456
1457 /*     FUNCTION : */
1458 /*     ---------- */
1459 /*     Normalization of constraints of derivatives, defined on DFUVIN */
1460 /*     on block DUVOUT. */
1461
1462 /*     KEYWORDS : */
1463 /*     ----------- */
1464 /*     ALL, AB_SPECIFI::VECTEUR&,DERIVEE&,NORMALISATION,&VECTEUR */
1465
1466 /*     INPUT ARGUMENTS : */
1467 /*     ------------------ */
1468 /*     DFUVIN: Limits of the block of definition by U and by V where 
1469 */
1470 /*             constraints CNTRIN are defined. */
1471 /*     NDIMEN: Dimension of the space. */
1472 /*     IORDRE: Order of constraint imposed at the extremities of the iso. */
1473 /*             (if Iso-U, it is necessary to calculate derivatives by V and vice */
1474 /*             versa). */
1475 /*             = 0, the extremities of the iso are calculated */
1476 /*             = 1, additionally the 1st derivative in the direction */
1477 /*                  of the iso is calculated */
1478 /*             = 2, additionally the 2nd derivative in the direction */
1479 /*                  of the iso is calculated */
1480 /*     CNTRIN: Contains, if IORDRE>=0, IORDRE+1 derivatives */
1481 /*             of order IORDRE of F(Uc,v) or of F(u,Vc), following the */
1482 /*             value of ISOFAV, RENORMALIZED by u and v in (-1,1). */
1483 /*     DUVOUT: Limits of the block of definition by U and by V where the */
1484 /*             constraints CNTOUT will be defined. */
1485 /*     ISOFAV: Isoparameter fixed for the discretization; */
1486 /*             = 1, discretization with fixed U=Uc and variable V. */
1487 /*             = 2, discretization with fixed V=Vc and variable U. */
1488 /*     IDERIV: Ordre de derivee transverse a l'iso fixee (Si Iso-U=Uc */
1489 /*             is fixed, the derivative of order IDERIV is discretized by U */
1490 /*             of F(Uc,v). The same if iso-V is fixed). */
1491 /*             Varies from (positioning) to 2 (2nd derivative). */
1492
1493 /*     OUTPUT ARGUMENTS : */
1494 /*     ------------------- */
1495 /*     CNTOUT: Contains, if IORDRE>=0, IORDRE+1 derivatives */
1496 /*             of order IORDRE of F(Uc,v) or of F(u,Vc), depending on the */
1497 /*             value of ISOFAV, RENORMALIZED for u and v in DUVOUT. */
1498
1499 /*     COMMONS USED   : */
1500 /*     ---------------- */
1501
1502 /*     REFERENCES CALLED   : */
1503 /*     --------------------- */
1504
1505 /*     DESCRIPTION/NOTES/LIMITATIONS : */
1506 /*     ------------------------------- */
1507 /*     CNTRIN can be an output/input  argument, */
1508 /*     so the call: */
1509
1510 /*      CALL MMA1NOC(DFUVIN,NDIMEN,IORDRE,CNTRIN,DUVOUT */
1511 /*     1           ,ISOFAV,IDERIV,CNTRIN) */
1512
1513 /*     is correct. */
1514 /* > */
1515 /* ********************************************************************** 
1516 */
1517
1518 /*   Name of the routine */
1519
1520
1521     /* Parameter adjustments */
1522     dfuvin -= 3;
1523     --cntout;
1524     --cntrin;
1525     duvout -= 3;
1526
1527     /* Function Body */
1528     ibb = AdvApp2Var_SysBase::mnfndeb_();
1529     if (ibb >= 3) {
1530         AdvApp2Var_SysBase::mgenmsg_("MMA1NOC", 7L);
1531     }
1532
1533 /* --------------- Determination of coefficients of normalization -------
1534  */
1535
1536     if (*isofav == 1) {
1537         d__1 = (dfuvin[4] - dfuvin[3]) / (duvout[4] - duvout[3]);
1538         rider = AdvApp2Var_MathBase::pow__di(&d__1, ideriv);
1539         d__1 = (dfuvin[6] - dfuvin[5]) / (duvout[6] - duvout[5]);
1540         riord = AdvApp2Var_MathBase::pow__di(&d__1, iordre);
1541
1542     } else {
1543         d__1 = (dfuvin[6] - dfuvin[5]) / (duvout[6] - duvout[5]);
1544         rider = AdvApp2Var_MathBase::pow__di(&d__1, ideriv);
1545         d__1 = (dfuvin[4] - dfuvin[3]) / (duvout[4] - duvout[3]);
1546         riord = AdvApp2Var_MathBase::pow__di(&d__1, iordre);
1547     }
1548
1549 /* ------------- Renormalization of the vector of constraint --------------- 
1550 */
1551
1552     bid = rider * riord;
1553     i__1 = *ndimen;
1554     for (nd = 1; nd <= i__1; ++nd) {
1555         cntout[nd] = bid * cntrin[nd];
1556 /* L100: */
1557     }
1558
1559 /* ------------------------------ The end ------------------------------- 
1560 */
1561
1562     if (ibb >= 3) {
1563         AdvApp2Var_SysBase::mgsomsg_("MMA1NOC", 7L);
1564     }
1565     return 0;
1566 } /* mma1noc_ */
1567
1568 //=======================================================================
1569 //function : mma1nop_
1570 //purpose  : 
1571 //=======================================================================
1572 int mma1nop_(integer *nbroot, 
1573              doublereal *rootlg, 
1574              doublereal *uvfonc, 
1575              integer *isofav, 
1576              doublereal *ttable, 
1577              integer *iercod)
1578
1579 {
1580   /* System generated locals */
1581   integer i__1;
1582
1583   /* Local variables */
1584   static doublereal alinu, blinu, alinv, blinv;
1585   static integer ii, ibb;
1586
1587 /* ***********************************************************************
1588  */
1589
1590 /*     FUNCTION : */
1591 /*     ---------- */
1592 /*     Normalization of parameters of an iso, starting from  */
1593 /*     parametric block and parameters on (-1,1). */
1594
1595 /*     KEYWORDS : */
1596 /*     ----------- */
1597 /*      TOUS,AB_SPECIFI :: ISO&,POINT&,NORMALISATION,&POINT,&ISO */
1598
1599 /*     INPUT ARGUMENTS : */
1600 /*     -------------------- */
1601 /*        NBROOT: Nb of points of discretisation INSIDE the iso */
1602 /*                defined on (-1,1). */
1603 /*        ROOTLG: Table of discretization parameters on )-1,1( */
1604 /*                of the iso. */
1605 /*        UVFONC: Block of definition of the iso */
1606 /*        ISOFAV: = 1, this is iso-u; =2, this is iso-v. */
1607
1608 /*     OUTPUT ARGUMENTS : */
1609 /*     --------------------- */
1610 /*        TTABLE: Table of parameters renormalized on UVFONC of the iso. 
1611 */
1612 /*        IERCOD: = 0, OK */
1613 /*                = 1, ISOFAV is out of allowed values. */
1614
1615 /* > */
1616 /* ********************************************************************** 
1617 */
1618 /*   Name of the routine */
1619
1620
1621     /* Parameter adjustments */
1622     --rootlg;
1623     uvfonc -= 3;
1624
1625     /* Function Body */
1626     ibb = AdvApp2Var_SysBase::mnfndeb_();
1627     if (ibb >= 3) {
1628         AdvApp2Var_SysBase::mgenmsg_("MMA1NOP", 7L);
1629     }
1630
1631     alinu = (uvfonc[4] - uvfonc[3]) / 2.;
1632     blinu = (uvfonc[4] + uvfonc[3]) / 2.;
1633     alinv = (uvfonc[6] - uvfonc[5]) / 2.;
1634     blinv = (uvfonc[6] + uvfonc[5]) / 2.;
1635
1636     if (*isofav == 1) {
1637         ttable[0] = uvfonc[5];
1638         i__1 = *nbroot;
1639         for (ii = 1; ii <= i__1; ++ii) {
1640             ttable[ii] = alinv * rootlg[ii] + blinv;
1641 /* L100: */
1642         }
1643         ttable[*nbroot + 1] = uvfonc[6];
1644     } else if (*isofav == 2) {
1645         ttable[0] = uvfonc[3];
1646         i__1 = *nbroot;
1647         for (ii = 1; ii <= i__1; ++ii) {
1648             ttable[ii] = alinu * rootlg[ii] + blinu;
1649 /* L200: */
1650         }
1651         ttable[*nbroot + 1] = uvfonc[4];
1652     } else {
1653         goto L9100;
1654     }
1655
1656     goto L9999;
1657
1658 /* ------------------------------ THE END ------------------------------- 
1659 */
1660
1661 L9100:
1662     *iercod = 1;
1663     goto L9999;
1664
1665 L9999:
1666     if (*iercod != 0) {
1667         AdvApp2Var_SysBase::maermsg_("MMA1NOP", iercod, 7L);
1668     }
1669     if (ibb >= 3) {
1670         AdvApp2Var_SysBase::mgsomsg_("MMA1NOP", 7L);
1671     }
1672
1673  return 0 ;
1674
1675 } /* mma1nop_ */
1676
1677 //=======================================================================
1678 //function : mma2ac1_
1679 //purpose  : 
1680 //=======================================================================
1681 int AdvApp2Var_ApproxF2var::mma2ac1_(integer const *ndimen, 
1682                                      integer const *mxujac, 
1683                                      integer const *mxvjac, 
1684                                      integer const *iordru, 
1685                                      integer const *iordrv, 
1686                                      doublereal const *contr1, 
1687                                      doublereal const * contr2, 
1688                                      doublereal const *contr3, 
1689                                      doublereal const *contr4, 
1690                                      doublereal const *uhermt, 
1691                                      doublereal const *vhermt, 
1692                                      doublereal *patjac)
1693
1694 {
1695   /* System generated locals */
1696   integer contr1_dim1, contr1_dim2, contr1_offset, contr2_dim1, contr2_dim2,
1697   contr2_offset, contr3_dim1, contr3_dim2, contr3_offset, 
1698   contr4_dim1, contr4_dim2, contr4_offset, uhermt_dim1, 
1699   uhermt_offset, vhermt_dim1, vhermt_offset, patjac_dim1, 
1700   patjac_dim2, patjac_offset, i__1, i__2, i__3, i__4, i__5;
1701
1702   /* Local variables */
1703   static logical ldbg;
1704   static integer ndgu, ndgv;
1705   static doublereal bidu1, bidu2, bidv1, bidv2;
1706   static integer ioru1, iorv1, ii, nd, jj, ku, kv;
1707   static doublereal cnt1, cnt2, cnt3, cnt4;
1708
1709 /* ********************************************************************** 
1710 */
1711
1712 /*     FUNCTION : */
1713 /*     ---------- */
1714 /*     Add polynoms of edge constraints. */
1715
1716 /*     KEYWORDS : */
1717 /*     ----------- */
1718 /*  TOUS,AB_SPECIFI::POINT&,CONTRAINTE&,ADDITION,&POLYNOME */
1719
1720 /*     INPUT ARGUMENTS  : */
1721 /*     ------------------ */
1722 /*   NDIMEN: Dimension of the space. */
1723 /*   MXUJAC: Max degree of the polynom of approximation by U. The  */
1724 /*           representation in the orthogonal base starts from degree */
1725 /*           0 to degree MXUJAC-2*(IORDRU+1). The polynomial base is the */
1726 /*           base of Jacobi of order -1 (Legendre), 0, 1 or 2 */
1727 /*   MXVJAC: Max degree of the polynom of approximation by V. The  */
1728 /*           representation in the orthogonal base starts from degree */
1729 /*           0 to degree MXUJAC-2*(IORDRU+1). The polynomial base is the */
1730 /*           base of Jacobi of order -1 (Legendre), 0, 1 or 2 */
1731 /*   IORDRU: Order of the base of Jacobi (-1,0,1 or 2) by U. Corresponds */
1732 /*           to the step of constraints: C0, C1 or C2. */
1733 /*   IORDRV: Order of the base of Jacobi (-1,0,1 or 2) by V. Corresponds */
1734 /*           to the step of constraints: C0, C1 or C2. */
1735 /*   CONTR1: Contains, if IORDRU and IORDRV>=0, the values at the */
1736 /*           extremities of F(U0,V0) and its derivatives. */
1737 /*   CONTR2: Contains, if IORDRU and IORDRV>=0, the values at the */
1738 /*           extremities of F(U1,V0) and its derivatives. */
1739 /*   CONTR3: Contains, if IORDRU and IORDRV>=0, the values at the */
1740 /*           extremities of F(U0,V1) and its derivatives. */
1741 /*   CONTR4: Contains, if IORDRU and IORDRV>=0, the values at the */
1742 /*           extremities of F(U1,V1) and its derivatives. */
1743 /*   UHERMT: Coeff. of Hermit polynoms of order IORDRU. */
1744 /*   VHERMT: Coeff. of Hermit polynoms of order IORDRV. */
1745 /*   PATJAC: Table of coefficients of the polynom P(u,v) of approximation */
1746 /*           of F(u,v) WITHOUT taking into account the constraints. */
1747
1748 /*     OUTPUT ARGUMENTS : */
1749 /*     ------------------- */
1750 /*   PATJAC: Table of coefficients of the polynom P(u,v) by approximation */
1751 /*           of F(u,v) WITH taking into account of constraints. */
1752 /* > */
1753 /* ********************************************************************** 
1754 */
1755 /*   Name of the routine */
1756
1757 /* --------------------------- Initialization -------------------------- 
1758 */
1759
1760     /* Parameter adjustments */
1761     patjac_dim1 = *mxujac + 1;
1762     patjac_dim2 = *mxvjac + 1;
1763     patjac_offset = patjac_dim1 * patjac_dim2;
1764     patjac -= patjac_offset;
1765     uhermt_dim1 = (*iordru << 1) + 2;
1766     uhermt_offset = uhermt_dim1;
1767     uhermt -= uhermt_offset;
1768     vhermt_dim1 = (*iordrv << 1) + 2;
1769     vhermt_offset = vhermt_dim1;
1770     vhermt -= vhermt_offset;
1771     contr4_dim1 = *ndimen;
1772     contr4_dim2 = *iordru + 2;
1773     contr4_offset = contr4_dim1 * (contr4_dim2 + 1) + 1;
1774     contr4 -= contr4_offset;
1775     contr3_dim1 = *ndimen;
1776     contr3_dim2 = *iordru + 2;
1777     contr3_offset = contr3_dim1 * (contr3_dim2 + 1) + 1;
1778     contr3 -= contr3_offset;
1779     contr2_dim1 = *ndimen;
1780     contr2_dim2 = *iordru + 2;
1781     contr2_offset = contr2_dim1 * (contr2_dim2 + 1) + 1;
1782     contr2 -= contr2_offset;
1783     contr1_dim1 = *ndimen;
1784     contr1_dim2 = *iordru + 2;
1785     contr1_offset = contr1_dim1 * (contr1_dim2 + 1) + 1;
1786     contr1 -= contr1_offset;
1787
1788     /* Function Body */
1789     ldbg = AdvApp2Var_SysBase::mnfndeb_() >= 3;
1790     if (ldbg) {
1791         AdvApp2Var_SysBase::mgenmsg_("MMA2AC1", 7L);
1792     }
1793
1794 /* ------------ SUBTRACTION OF ANGULAR CONSTRAINTS ------------------- 
1795 */
1796
1797     ioru1 = *iordru + 1;
1798     iorv1 = *iordrv + 1;
1799     ndgu = (*iordru << 1) + 1;
1800     ndgv = (*iordrv << 1) + 1;
1801
1802     i__1 = iorv1;
1803     for (jj = 1; jj <= i__1; ++jj) {
1804         i__2 = ioru1;
1805         for (ii = 1; ii <= i__2; ++ii) {
1806             i__3 = *ndimen;
1807             for (nd = 1; nd <= i__3; ++nd) {
1808                 cnt1 = contr1[nd + (ii + jj * contr1_dim2) * contr1_dim1];
1809                 cnt2 = contr2[nd + (ii + jj * contr2_dim2) * contr2_dim1];
1810                 cnt3 = contr3[nd + (ii + jj * contr3_dim2) * contr3_dim1];
1811                 cnt4 = contr4[nd + (ii + jj * contr4_dim2) * contr4_dim1];
1812                 i__4 = ndgv;
1813                 for (kv = 0; kv <= i__4; ++kv) {
1814                     bidv1 = vhermt[kv + ((jj << 1) - 1) * vhermt_dim1];
1815                     bidv2 = vhermt[kv + (jj << 1) * vhermt_dim1];
1816                     i__5 = ndgu;
1817                     for (ku = 0; ku <= i__5; ++ku) {
1818                         bidu1 = uhermt[ku + ((ii << 1) - 1) * uhermt_dim1];
1819                         bidu2 = uhermt[ku + (ii << 1) * uhermt_dim1];
1820                         patjac[ku + (kv + nd * patjac_dim2) * patjac_dim1] = 
1821                                 patjac[ku + (kv + nd * patjac_dim2) * 
1822                                 patjac_dim1] - bidu1 * bidv1 * cnt1 - bidu2 * 
1823                                 bidv1 * cnt2 - bidu1 * bidv2 * cnt3 - bidu2 * 
1824                                 bidv2 * cnt4;
1825 /* L500: */
1826                     }
1827 /* L400: */
1828                 }
1829 /* L300: */
1830             }
1831 /* L200: */
1832         }
1833 /* L100: */
1834     }
1835
1836 /* ------------------------------ The end ------------------------------- 
1837 */
1838
1839     if (ldbg) {
1840         AdvApp2Var_SysBase::mgsomsg_("MMA2AC1", 7L);
1841     }
1842     return 0;
1843 } /* mma2ac1_ */
1844
1845 //=======================================================================
1846 //function : mma2ac2_
1847 //purpose  : 
1848 //=======================================================================
1849 int AdvApp2Var_ApproxF2var::mma2ac2_(const integer *ndimen, 
1850                                      const integer *mxujac, 
1851                                      const integer *mxvjac, 
1852                                      const integer *iordrv, 
1853                                      const integer *nclimu, 
1854                                      const integer *ncfiv1, 
1855                                      const doublereal *crbiv1, 
1856                                      const integer *ncfiv2, 
1857                                      const doublereal *crbiv2, 
1858                                      const doublereal *vhermt, 
1859                                      doublereal *patjac)
1860
1861 {
1862   /* System generated locals */
1863   integer crbiv1_dim1, crbiv1_dim2, crbiv1_offset, crbiv2_dim1, crbiv2_dim2,
1864   crbiv2_offset, patjac_dim1, patjac_dim2, patjac_offset, 
1865   vhermt_dim1, vhermt_offset, i__1, i__2, i__3, i__4;
1866
1867   /* Local variables */
1868   static logical ldbg;
1869   static integer ndgv1, ndgv2, ii, jj, nd, kk;
1870   static doublereal bid1, bid2;
1871
1872 /* ********************************************************************** 
1873 */
1874
1875 /*     FUNCTION : */
1876 /*     ---------- */
1877 /*     Add polynoms of constraints */
1878
1879 /*     KEYWORDS : */
1880 /*     ----------- */
1881 /*     FUNCTION,APPROXIMATION,COEFFICIENT,POLYNOM */
1882
1883 /*     INPUT ARGUMENTS : */
1884 /*     ------------------ */
1885 /*   NDIMEN: Dimension of the space. */
1886 /*   MXUJAC: Max degree of the polynom of approximation by U. The  */
1887 /*           representation in the orthogonal base starts from degree */
1888 /*           0 to degree MXUJAC-2*(IORDRU+1). The polynomial base is the */
1889 /*           base of Jacobi of order -1 (Legendre), 0, 1 or 2 */
1890 /*   MXVJAC: Max degree of the polynom of approximation by V. The  */
1891 /*           representation in the orthogonal base starts from degree */
1892 /*           0 to degree MXUJAC-2*(IORDRU+1). The polynomial base is the */
1893 /*           base of Jacobi of order -1 (Legendre), 0, 1 or 2 */
1894 /*   IORDRV: Order of the base of Jacobi (-1,0,1 or 2) by V. Corresponds */
1895 /*           to the step of constraints: C0, C1 or C2. */
1896 /*   NCLIMU  LIMIT nb of coeff by u of the solution P(u,v) 
1897 *    NCFIV1: Nb of Coeff. of curves stored in CRBIV1. */
1898 /*   CRBIV1: Table of coeffs of the approximation of iso-V0 and its */
1899 /*           derivatives till order IORDRV. */
1900 /*   NCFIV2: Nb of Coeff. of curves stored in CRBIV2. */
1901 /*   CRBIV2: Table of coeffs of approximation of iso-V1 and its */
1902 /*           derivatives till order IORDRV. */
1903 /*   VHERMT: Coeff. of Hermit polynoms of order IORDRV. */
1904 /*   PATJAC: Table of coefficients of the polynom P(u,v) of approximation */
1905 /*           of F(u,v) WITHOUT taking into account the constraints. */
1906
1907 /*     OUTPUT ARGUMENTS : */
1908 /*     ------------------- */
1909 /*   PATJAC: Table of coefficients of the polynom P(u,v) by approximation */
1910 /*           of F(u,v) WITH taking into account of constraints. */
1911 /* > *//*
1912
1913
1914 /* > */
1915 /* ********************************************************************** 
1916 */
1917 /*   Name of the routine */
1918
1919 /* --------------------------- Initialisations -------------------------- 
1920 */
1921
1922     /* Parameter adjustments */
1923     patjac_dim1 = *mxujac + 1;
1924     patjac_dim2 = *mxvjac + 1;
1925     patjac_offset = patjac_dim1 * patjac_dim2;
1926     patjac -= patjac_offset;
1927     vhermt_dim1 = (*iordrv << 1) + 2;
1928     vhermt_offset = vhermt_dim1;
1929     vhermt -= vhermt_offset;
1930     --ncfiv2;
1931     --ncfiv1;
1932     crbiv2_dim1 = *nclimu;
1933     crbiv2_dim2 = *ndimen;
1934     crbiv2_offset = crbiv2_dim1 * (crbiv2_dim2 + 1);
1935     crbiv2 -= crbiv2_offset;
1936     crbiv1_dim1 = *nclimu;
1937     crbiv1_dim2 = *ndimen;
1938     crbiv1_offset = crbiv1_dim1 * (crbiv1_dim2 + 1);
1939     crbiv1 -= crbiv1_offset;
1940
1941     /* Function Body */
1942     ldbg = AdvApp2Var_SysBase::mnfndeb_() >= 3;
1943     if (ldbg) {
1944         AdvApp2Var_SysBase::mgenmsg_("MMA2AC2", 7L);
1945     }
1946
1947 /* ------------ ADDING of coeff by u of curves, by v of Hermit -------- 
1948 */
1949
1950     i__1 = *iordrv + 1;
1951     for (ii = 1; ii <= i__1; ++ii) {
1952         ndgv1 = ncfiv1[ii] - 1;
1953         ndgv2 = ncfiv2[ii] - 1;
1954         i__2 = *ndimen;
1955         for (nd = 1; nd <= i__2; ++nd) {
1956             i__3 = (*iordrv << 1) + 1;
1957             for (jj = 0; jj <= i__3; ++jj) {
1958                 bid1 = vhermt[jj + ((ii << 1) - 1) * vhermt_dim1];
1959                 i__4 = ndgv1;
1960                 for (kk = 0; kk <= i__4; ++kk) {
1961                     patjac[kk + (jj + nd * patjac_dim2) * patjac_dim1] += 
1962                             bid1 * crbiv1[kk + (nd + ii * crbiv1_dim2) * 
1963                             crbiv1_dim1];
1964 /* L400: */
1965                 }
1966                 bid2 = vhermt[jj + (ii << 1) * vhermt_dim1];
1967                 i__4 = ndgv2;
1968                 for (kk = 0; kk <= i__4; ++kk) {
1969                     patjac[kk + (jj + nd * patjac_dim2) * patjac_dim1] += 
1970                             bid2 * crbiv2[kk + (nd + ii * crbiv2_dim2) * 
1971                             crbiv2_dim1];
1972 /* L500: */
1973                 }
1974 /* L300: */
1975             }
1976 /* L200: */
1977         }
1978 /* L100: */
1979     }
1980
1981 /* ------------------------------ The end ------------------------------- 
1982 */
1983
1984     if (ldbg) {
1985         AdvApp2Var_SysBase::mgsomsg_("MMA2AC2", 7L);
1986     }
1987     return 0;
1988 } /* mma2ac2_ */
1989
1990
1991 //=======================================================================
1992 //function : mma2ac3_
1993 //purpose  : 
1994 //=======================================================================
1995 int AdvApp2Var_ApproxF2var::mma2ac3_(const integer *ndimen, 
1996                                      const integer *mxujac, 
1997                                      const integer *mxvjac, 
1998                                      const integer *iordru, 
1999                                      const integer *nclimv, 
2000                                      const integer *ncfiu1, 
2001                                      const doublereal * crbiu1, 
2002                                      const integer *ncfiu2, 
2003                                      const doublereal *crbiu2, 
2004                                      const doublereal *uhermt, 
2005                                      doublereal *patjac)
2006
2007 {
2008   /* System generated locals */
2009   integer crbiu1_dim1, crbiu1_dim2, crbiu1_offset, crbiu2_dim1, crbiu2_dim2,
2010   crbiu2_offset, patjac_dim1, patjac_dim2, patjac_offset, 
2011   uhermt_dim1, uhermt_offset, i__1, i__2, i__3, i__4;
2012
2013   /* Local variables */
2014   static logical ldbg;
2015   static integer ndgu1, ndgu2, ii, jj, nd, kk;
2016   static doublereal bid1, bid2;
2017
2018 /* ********************************************************************** 
2019 */
2020
2021 /*     FUNCTION : */
2022 /*     ---------- */
2023 /*     Ajout des polynomes de contraintes */
2024
2025 /*     KEYWORDS : */
2026 /*     ----------- */
2027 /*     FONCTION,APPROXIMATION,COEFFICIENT,POLYNOME */
2028
2029 /*     INPUT ARGUMENTS : */
2030 /*     ------------------ */
2031 /*   NDIMEN: Dimension of the space. */
2032 /*   MXUJAC: Max degree of the polynom of approximation by U. The  */
2033 /*           representation in the orthogonal base starts from degree */
2034 /*           0 to degree MXUJAC-2*(IORDRU+1). The polynomial base is the */
2035 /*           base of Jacobi of order -1 (Legendre), 0, 1 or 2 */
2036 /*   MXVJAC: Max degree of the polynom of approximation by V. The  */
2037 /*           representation in the orthogonal base starts from degree */
2038 /*           0 to degree MXUJAC-2*(IORDRU+1). The polynomial base is the */
2039 /*           base of Jacobi of order -1 (Legendre), 0, 1 or 2 */
2040 /*   IORDRU: Order of the base of Jacobi (-1,0,1 or 2) by U. Corresponds */
2041 /*           to the step of constraints: C0, C1 or C2. */
2042 /*   NCLIMV  LIMIT nb of coeff by v of the solution P(u,v) 
2043 *    NCFIU1: Nb of Coeff. of curves stored in CRBIU1. */
2044 /*   CRBIU1: Table of coeffs of the approximation of iso-U0 and its */
2045 /*           derivatives till order IORDRU. */
2046 /*   NCFIU2: Nb of Coeff. of curves stored in CRBIU2. */
2047 /*   CRBIU2: Table of coeffs of approximation of iso-U1 and its */
2048 /*           derivatives till order IORDRU */
2049 /*   UHERMT: Coeff. of Hermit polynoms of order IORDRU. */
2050 /*   PATJAC: Table of coefficients of the polynom P(u,v) of approximation */
2051 /*           of F(u,v) WITHOUT taking into account the constraints. */
2052
2053 /*     OUTPUT ARGUMENTS : */
2054 /*     ------------------- */
2055 /*   PATJAC: Table of coefficients of the polynom P(u,v) by approximation */
2056 /*           of F(u,v) WITH taking into account of constraints. */
2057
2058
2059 /* > */
2060 /* ********************************************************************** 
2061 */
2062 /*   The name of the routine */
2063
2064 /* --------------------------- Initializations -------------------------- 
2065 */
2066
2067     /* Parameter adjustments */
2068     patjac_dim1 = *mxujac + 1;
2069     patjac_dim2 = *mxvjac + 1;
2070     patjac_offset = patjac_dim1 * patjac_dim2;
2071     patjac -= patjac_offset;
2072     uhermt_dim1 = (*iordru << 1) + 2;
2073     uhermt_offset = uhermt_dim1;
2074     uhermt -= uhermt_offset;
2075     --ncfiu2;
2076     --ncfiu1;
2077     crbiu2_dim1 = *nclimv;
2078     crbiu2_dim2 = *ndimen;
2079     crbiu2_offset = crbiu2_dim1 * (crbiu2_dim2 + 1);
2080     crbiu2 -= crbiu2_offset;
2081     crbiu1_dim1 = *nclimv;
2082     crbiu1_dim2 = *ndimen;
2083     crbiu1_offset = crbiu1_dim1 * (crbiu1_dim2 + 1);
2084     crbiu1 -= crbiu1_offset;
2085
2086     /* Function Body */
2087     ldbg = AdvApp2Var_SysBase::mnfndeb_() >= 3;
2088     if (ldbg) {
2089         AdvApp2Var_SysBase::mgenmsg_("MMA2AC3", 7L);
2090     }
2091
2092 /* ------------ ADDING of coeff by u of curves, by v of Hermit -------- 
2093 */
2094
2095     i__1 = *iordru + 1;
2096     for (ii = 1; ii <= i__1; ++ii) {
2097         ndgu1 = ncfiu1[ii] - 1;
2098         ndgu2 = ncfiu2[ii] - 1;
2099         i__2 = *ndimen;
2100         for (nd = 1; nd <= i__2; ++nd) {
2101             i__3 = ndgu1;
2102             for (jj = 0; jj <= i__3; ++jj) {
2103                 bid1 = crbiu1[jj + (nd + ii * crbiu1_dim2) * crbiu1_dim1];
2104                 i__4 = (*iordru << 1) + 1;
2105                 for (kk = 0; kk <= i__4; ++kk) {
2106                     patjac[kk + (jj + nd * patjac_dim2) * patjac_dim1] += 
2107                             bid1 * uhermt[kk + ((ii << 1) - 1) * uhermt_dim1];
2108 /* L400: */
2109                 }
2110 /* L300: */
2111             }
2112             i__3 = ndgu2;
2113             for (jj = 0; jj <= i__3; ++jj) {
2114                 bid2 = crbiu2[jj + (nd + ii * crbiu2_dim2) * crbiu2_dim1];
2115                 i__4 = (*iordru << 1) + 1;
2116                 for (kk = 0; kk <= i__4; ++kk) {
2117                     patjac[kk + (jj + nd * patjac_dim2) * patjac_dim1] += 
2118                             bid2 * uhermt[kk + (ii << 1) * uhermt_dim1];
2119 /* L600: */
2120                 }
2121 /* L500: */
2122             }
2123
2124 /* L200: */
2125         }
2126 /* L100: */
2127     }
2128
2129 /* ------------------------------ The end ------------------------------- 
2130 */
2131
2132     if (ldbg) {
2133         AdvApp2Var_SysBase::mgsomsg_("MMA2AC3", 7L);
2134     }
2135     return 0;
2136 } /* mma2ac3_ */
2137
2138 //=======================================================================
2139 //function : mma2can_
2140 //purpose  : 
2141 //=======================================================================
2142 int AdvApp2Var_ApproxF2var::mma2can_(const integer *ncfmxu, 
2143                                      const integer *ncfmxv,
2144                                      const integer *ndimen, 
2145                                      const integer *iordru, 
2146                                      const integer *iordrv, 
2147                                      const integer *ncoefu, 
2148                                      const integer *ncoefv, 
2149                                      const doublereal *patjac, 
2150                                      doublereal *pataux, 
2151                                      doublereal *patcan, 
2152                                      integer *iercod)
2153
2154 {
2155   /* System generated locals */
2156   integer patjac_dim1, patjac_dim2, patjac_offset, patcan_dim1, patcan_dim2,
2157   patcan_offset, i__1, i__2;
2158
2159   /* Local variables */
2160   static logical ldbg;
2161   static integer ilon1, ilon2, ii, nd;
2162
2163 /* ********************************************************************** 
2164 */
2165
2166 /*     FUNCTION : */
2167 /*     ---------- */
2168 /*     Change of Jacobi base to canonical (-1,1) and writing in a greater */
2169 /*     table. */
2170
2171 /*     KEYWORDS : */
2172 /*     ----------- */
2173 /*     ALL,AB_SPECIFI,CARREAU&,CONVERSION,JACOBI,CANNONIQUE,&CARREAU */
2174
2175 /*     INPUT ARGUMENTS : */
2176 /*     -------------------- */
2177 /*     NCFMXU: Dimension by U of resulting table PATCAN */
2178 /*     NCFMXV: Dimension by V of resulting table PATCAN */
2179 /*     NDIMEN: Dimension of the workspace. */
2180 /*     IORDRU: Order of constraint by U */
2181 /*     IORDRV: Order of constraint by V. */
2182 /*     NCOEFU: Nb of coeff by U of square PATJAC */
2183 /*     NCOEFV: Nb of coeff by V of square PATJAC */
2184 /*     PATJAC: Square in the base of Jacobi of order IORDRU by U and */
2185 /*             IORDRV by V. */
2186
2187 /*     OUTPUT ARGUMENTS : */
2188 /*     --------------------- */
2189 /*     PATAUX: Auxiliary Table. */
2190 /*     PATCAN: Table of coefficients in the canonic base. */
2191 /*     IERCOD: Error code. */
2192 /*             = 0, everything goes well, and all things are equal. */
2193 /*             = 1, the program refuses to process with incorrect input arguments */
2194
2195
2196 /*     COMMONS USED : */
2197 /*     ------------------ */
2198
2199 /*     REFERENCES CALLED : */
2200 /*     --------------------- */
2201
2202 /*     DESCRIPTION/NOTES/LIMITATIONS : */
2203 /*     ----------------------------------- */
2204 /* > */
2205 /* ********************************************************************** 
2206 */
2207
2208
2209     /* Parameter adjustments */
2210     patcan_dim1 = *ncfmxu;
2211     patcan_dim2 = *ncfmxv;
2212     patcan_offset = patcan_dim1 * (patcan_dim2 + 1) + 1;
2213     patcan -= patcan_offset;
2214     --pataux;
2215     patjac_dim1 = *ncoefu;
2216     patjac_dim2 = *ncoefv;
2217     patjac_offset = patjac_dim1 * (patjac_dim2 + 1) + 1;
2218     patjac -= patjac_offset;
2219
2220     /* Function Body */
2221     ldbg = AdvApp2Var_SysBase::mnfndeb_() >= 2;
2222     if (ldbg) {
2223         AdvApp2Var_SysBase::mgenmsg_("MMA2CAN", 7L);
2224     }
2225     *iercod = 0;
2226
2227     if (*iordru < -1 || *iordru > 2) {
2228         goto L9100;
2229     }
2230     if (*iordrv < -1 || *iordrv > 2) {
2231         goto L9100;
2232     }
2233     if (*ncoefu > *ncfmxu || *ncoefv > *ncfmxv) {
2234         goto L9100;
2235     }
2236
2237 /* --> Pass to canonic base (-1,1) */
2238     mmjacpt_(ndimen, ncoefu, ncoefv, iordru, iordrv, &patjac[patjac_offset], &
2239             pataux[1], &patcan[patcan_offset]);
2240
2241 /* --> Write all in a greater table */
2242     AdvApp2Var_MathBase::mmfmca8_(ncoefu, 
2243              ncoefv, 
2244              ndimen, 
2245              ncfmxu, 
2246              ncfmxv, 
2247              ndimen, 
2248              &patcan[patcan_offset], 
2249              &patcan[patcan_offset]);
2250
2251 /* --> Complete with zeros the resulting table. */
2252     ilon1 = *ncfmxu - *ncoefu;
2253     ilon2 = *ncfmxu * (*ncfmxv - *ncoefv);
2254     i__1 = *ndimen;
2255     for (nd = 1; nd <= i__1; ++nd) {
2256         if (ilon1 > 0) {
2257             i__2 = *ncoefv;
2258             for (ii = 1; ii <= i__2; ++ii) {
2259                 AdvApp2Var_SysBase::mvriraz_(&ilon1, 
2260                          &patcan[*ncoefu + 1 + (ii + nd * patcan_dim2) * patcan_dim1]);
2261 /* L110: */
2262             }
2263         }
2264         if (ilon2 > 0) {
2265             AdvApp2Var_SysBase::mvriraz_(&ilon2, 
2266                      &patcan[(*ncoefv + 1 + nd * patcan_dim2) * patcan_dim1 + 1]);
2267         }
2268 /* L100: */
2269     }
2270
2271     goto L9999;
2272
2273 /* ---------------------- 
2274 */
2275
2276 L9100:
2277     *iercod = 1;
2278     goto L9999;
2279
2280 L9999:
2281     AdvApp2Var_SysBase::maermsg_("MMA2CAN", iercod, 7L);
2282     if (ldbg) {
2283         AdvApp2Var_SysBase::mgsomsg_("MMA2CAN", 7L);
2284     }
2285  return 0 ;
2286 } /* mma2can_ */
2287
2288 //=======================================================================
2289 //function : mma2cd1_
2290 //purpose  : 
2291 //=======================================================================
2292 int mma2cd1_(integer *ndimen, 
2293              integer *nbpntu, 
2294              doublereal *urootl, 
2295              integer *nbpntv, 
2296              doublereal *vrootl, 
2297              integer *iordru, 
2298              integer *iordrv, 
2299              doublereal *contr1, 
2300              doublereal *contr2, 
2301              doublereal *contr3, 
2302              doublereal *contr4, 
2303              doublereal *fpntbu, 
2304              doublereal *fpntbv, 
2305              doublereal *uhermt, 
2306              doublereal *vhermt, 
2307              doublereal *sosotb, 
2308              doublereal *soditb, 
2309              doublereal *disotb, 
2310              doublereal *diditb)
2311
2312 {
2313   static integer c__1 = 1;
2314
2315 /* System generated locals */
2316     integer contr1_dim1, contr1_dim2, contr1_offset, contr2_dim1, contr2_dim2,
2317              contr2_offset, contr3_dim1, contr3_dim2, contr3_offset, 
2318             contr4_dim1, contr4_dim2, contr4_offset, uhermt_dim1, 
2319             uhermt_offset, vhermt_dim1, vhermt_offset, fpntbu_dim1, 
2320             fpntbu_offset, fpntbv_dim1, fpntbv_offset, sosotb_dim1, 
2321             sosotb_dim2, sosotb_offset, diditb_dim1, diditb_dim2, 
2322             diditb_offset, soditb_dim1, soditb_dim2, soditb_offset, 
2323             disotb_dim1, disotb_dim2, disotb_offset, i__1, i__2, i__3, i__4, 
2324             i__5;
2325
2326     /* Local variables */
2327     static integer ncfhu, ncfhv, nuroo, nvroo, nd, ii, jj, kk, ll, ibb, kkm, 
2328             llm, kkp, llp;
2329     static doublereal bid1, bid2, bid3, bid4;
2330     static doublereal diu1, diu2, div1, div2, sou1, sou2, sov1, sov2;
2331
2332 /* ********************************************************************** 
2333 */
2334
2335 /*     FUNCTION : */
2336 /*     ---------- */
2337 /*     Discretisation on the parameters of polynoms of interpolation */
2338 /*     of constraints at the corners of order IORDRE. */
2339
2340 /*     KEYWORDS : */
2341 /*     ----------- */
2342 /*     TOUS, AB_SPECIFI::CONTRAINTE&, DISCRETISATION, &POINT */
2343
2344 /*     INPUT ARGUMENTS : */
2345 /*     ------------------ */
2346 /*     NDIMEN: Dimension of the space. */
2347 /*     NBPNTU: Nb of INTERNAL parameters of discretisation by U. */
2348 /*             This is also the nb of root of Legendre polynom where discretization is done. */
2349 /*     UROOTL: Table of parameters of discretisation ON (-1,1) by U. 
2350 */
2351 /*     NBPNTV: Nb of INTERNAL  parameters of discretisation by V. */
2352 /*             This is also the nb of root of Legendre polynom where discretization is done. */
2353 /*     VROOTL: Table of discretization parameters on (-1,1) by V. 
2354 /*     IORDRU: Order of constraint imposed at the extremities of iso-V */
2355 /*             = 0, calculate the extremities of iso-V */
2356 /*             = 1, calculate, additionally, the 1st derivative in the direction of iso-V */
2357 /*             = 2, calculate, additionally, the 2nd derivative in the direction of iso-V */
2358 /*     IORDRV: Order of constraint imposed at the extremities of iso-U */
2359 /*             = 0, calculate the extremities of iso-U */
2360 /*             = 1, calculate, additionally, the 1st derivative in the direction of iso-U */
2361 /*             = 2, calculate, additionally, the 2nd derivative in the direction of iso-U */
2362 /*   CONTR1: Contains, if IORDRU and IORDRV>=0, the values at the */
2363 /*           extremities of F(U0,V0) and its derivatives. */
2364 /*   CONTR2: Contains, if IORDRU and IORDRV>=0, the values at the */
2365 /*           extremities of F(U1,V0) and its derivatives. */
2366 /*   CONTR3: Contains, if IORDRU and IORDRV>=0, the values at the */
2367 /*           extremities of F(U0,V1) and its derivatives. */
2368 /*   CONTR4: Contains, if IORDRU and IORDRV>=0, the values at the */
2369 /*           extremities of F(U1,V1) and its derivatives. */
2370 /*     SOSOTB: Preinitialized table (input/output argument). */
2371 /*     DISOTB: Preinitialized table (input/output argument). */
2372 /*     SODITB: Preinitialized table (input/output argument). */
2373 /*     DIDITB: Preinitialized table (input/output argument) */
2374
2375 /*     OUTPUT ARGUMENTS : */
2376 /*     ------------------- */
2377 /*     FPNTBU: Auxiliary table. */
2378 /*     FPNTBV: Auxiliary table. */
2379 /*     UHERMT: Table of 2*(IORDRU+1) coeff. of 2*(IORDRU+1) polynoms of Hermite. */
2380 /*     VHERMT: Table of 2*(IORDRV+1) coeff. of 2*(IORDRV+1) polynoms of Hermite. */
2381 /*   SOSOTB: Table where the terms of constraints are added */
2382 /*           C(ui,vj) + C(ui,-vj) + C(-ui,vj) + C(-ui,-vj) */
2383 /*           with ui and vj positive roots of the Legendre polynom */
2384 /*           of degree NBPNTU and NBPNTV respectively. */
2385 /*   DISOTB: Table where the terms of constraints are added */
2386 /*           C(ui,vj) + C(ui,-vj) - C(-ui,vj) - C(-ui,-vj) */
2387 /*           with ui and vj positive roots of the polynom of Legendre */
2388 /*           of degree NBPNTU and NBPNTV respectively. */
2389 /*   SODITB: Table where the terms of constraints are added */
2390 /*           C(ui,vj) - C(ui,-vj) + C(-ui,vj) - C(-ui,-vj) */
2391 /*           with ui and vj positive roots of the polynom of Legendre */
2392 /*           of degree NBPNTU and NBPNTV respectively. */
2393 /*   DIDITB: Table where the terms of constraints are added */
2394 /*           C(ui,vj) - C(ui,-vj) - C(-ui,vj) + C(-ui,-vj) */
2395 /*           with ui and vj positive roots of the polynom of Legendre */
2396 /*           of degree NBPNTU and NBPNTV respectively. */
2397
2398 /*     COMMONS USED   : */
2399 /*     ---------------- */
2400
2401 /*     REFERENCES CALLED   : */
2402 /*     ----------------------- */
2403
2404 /*     DESCRIPTION/NOTES/LIMITATIONS : */
2405 /*     ----------------------------------- */
2406
2407 /* > */
2408 /* ********************************************************************** 
2409 */
2410
2411 /*   Name of the routine */
2412
2413
2414     /* Parameter adjustments */
2415     --urootl;
2416     diditb_dim1 = *nbpntu / 2 + 1;
2417     diditb_dim2 = *nbpntv / 2 + 1;
2418     diditb_offset = diditb_dim1 * diditb_dim2;
2419     diditb -= diditb_offset;
2420     disotb_dim1 = *nbpntu / 2;
2421     disotb_dim2 = *nbpntv / 2;
2422     disotb_offset = disotb_dim1 * (disotb_dim2 + 1) + 1;
2423     disotb -= disotb_offset;
2424     soditb_dim1 = *nbpntu / 2;
2425     soditb_dim2 = *nbpntv / 2;
2426     soditb_offset = soditb_dim1 * (soditb_dim2 + 1) + 1;
2427     soditb -= soditb_offset;
2428     sosotb_dim1 = *nbpntu / 2 + 1;
2429     sosotb_dim2 = *nbpntv / 2 + 1;
2430     sosotb_offset = sosotb_dim1 * sosotb_dim2;
2431     sosotb -= sosotb_offset;
2432     --vrootl;
2433     uhermt_dim1 = (*iordru << 1) + 2;
2434     uhermt_offset = uhermt_dim1;
2435     uhermt -= uhermt_offset;
2436     fpntbu_dim1 = *nbpntu;
2437     fpntbu_offset = fpntbu_dim1 + 1;
2438     fpntbu -= fpntbu_offset;
2439     vhermt_dim1 = (*iordrv << 1) + 2;
2440     vhermt_offset = vhermt_dim1;
2441     vhermt -= vhermt_offset;
2442     fpntbv_dim1 = *nbpntv;
2443     fpntbv_offset = fpntbv_dim1 + 1;
2444     fpntbv -= fpntbv_offset;
2445     contr4_dim1 = *ndimen;
2446     contr4_dim2 = *iordru + 2;
2447     contr4_offset = contr4_dim1 * (contr4_dim2 + 1) + 1;
2448     contr4 -= contr4_offset;
2449     contr3_dim1 = *ndimen;
2450     contr3_dim2 = *iordru + 2;
2451     contr3_offset = contr3_dim1 * (contr3_dim2 + 1) + 1;
2452     contr3 -= contr3_offset;
2453     contr2_dim1 = *ndimen;
2454     contr2_dim2 = *iordru + 2;
2455     contr2_offset = contr2_dim1 * (contr2_dim2 + 1) + 1;
2456     contr2 -= contr2_offset;
2457     contr1_dim1 = *ndimen;
2458     contr1_dim2 = *iordru + 2;
2459     contr1_offset = contr1_dim1 * (contr1_dim2 + 1) + 1;
2460     contr1 -= contr1_offset;
2461
2462     /* Function Body */
2463     ibb = AdvApp2Var_SysBase::mnfndeb_();
2464     if (ibb >= 3) {
2465         AdvApp2Var_SysBase::mgenmsg_("MMA2CD1", 7L);
2466     }
2467
2468 /* ------------------- Discretisation of Hermite polynoms ----------- 
2469 */
2470
2471     ncfhu = (*iordru + 1) << 1;
2472     i__1 = ncfhu;
2473     for (ii = 1; ii <= i__1; ++ii) {
2474         i__2 = *nbpntu;
2475         for (ll = 1; ll <= i__2; ++ll) {
2476             AdvApp2Var_MathBase::mmmpocur_(&ncfhu, &c__1, &ncfhu, &uhermt[ii * uhermt_dim1], &
2477                     urootl[ll], &fpntbu[ll + ii * fpntbu_dim1]);
2478 /* L20: */
2479         }
2480 /* L10: */
2481     }
2482     ncfhv = (*iordrv + 1) << 1;
2483     i__1 = ncfhv;
2484     for (jj = 1; jj <= i__1; ++jj) {
2485         i__2 = *nbpntv;
2486         for (kk = 1; kk <= i__2; ++kk) {
2487             AdvApp2Var_MathBase::mmmpocur_(&ncfhv, &c__1, &ncfhv, &vhermt[jj * vhermt_dim1], &
2488                     vrootl[kk], &fpntbv[kk + jj * fpntbv_dim1]);
2489 /* L40: */
2490         }
2491 /* L30: */
2492     }
2493
2494 /* ---- The discretizations of polynoms of constraints are subtracted ---- 
2495 */
2496
2497     nuroo = *nbpntu / 2;
2498     nvroo = *nbpntv / 2;
2499     i__1 = *ndimen;
2500     for (nd = 1; nd <= i__1; ++nd) {
2501
2502         i__2 = *iordrv + 1;
2503         for (jj = 1; jj <= i__2; ++jj) {
2504             i__3 = *iordru + 1;
2505             for (ii = 1; ii <= i__3; ++ii) {
2506                 bid1 = contr1[nd + (ii + jj * contr1_dim2) * contr1_dim1];
2507                 bid2 = contr2[nd + (ii + jj * contr2_dim2) * contr2_dim1];
2508                 bid3 = contr3[nd + (ii + jj * contr3_dim2) * contr3_dim1];
2509                 bid4 = contr4[nd + (ii + jj * contr4_dim2) * contr4_dim1];
2510
2511                 i__4 = nvroo;
2512                 for (kk = 1; kk <= i__4; ++kk) {
2513                     kkp = (*nbpntv + 1) / 2 + kk;
2514                     kkm = nvroo - kk + 1;
2515                     sov1 = fpntbv[kkp + ((jj << 1) - 1) * fpntbv_dim1] + 
2516                             fpntbv[kkm + ((jj << 1) - 1) * fpntbv_dim1];
2517                     div1 = fpntbv[kkp + ((jj << 1) - 1) * fpntbv_dim1] - 
2518                             fpntbv[kkm + ((jj << 1) - 1) * fpntbv_dim1];
2519                     sov2 = fpntbv[kkp + (jj << 1) * fpntbv_dim1] + fpntbv[kkm 
2520                             + (jj << 1) * fpntbv_dim1];
2521                     div2 = fpntbv[kkp + (jj << 1) * fpntbv_dim1] - fpntbv[kkm 
2522                             + (jj << 1) * fpntbv_dim1];
2523                     i__5 = nuroo;
2524                     for (ll = 1; ll <= i__5; ++ll) {
2525                         llp = (*nbpntu + 1) / 2 + ll;
2526                         llm = nuroo - ll + 1;
2527                         sou1 = fpntbu[llp + ((ii << 1) - 1) * fpntbu_dim1] + 
2528                                 fpntbu[llm + ((ii << 1) - 1) * fpntbu_dim1];
2529                         diu1 = fpntbu[llp + ((ii << 1) - 1) * fpntbu_dim1] - 
2530                                 fpntbu[llm + ((ii << 1) - 1) * fpntbu_dim1];
2531                         sou2 = fpntbu[llp + (ii << 1) * fpntbu_dim1] + fpntbu[
2532                                 llm + (ii << 1) * fpntbu_dim1];
2533                         diu2 = fpntbu[llp + (ii << 1) * fpntbu_dim1] - fpntbu[
2534                                 llm + (ii << 1) * fpntbu_dim1];
2535                         sosotb[ll + (kk + nd * sosotb_dim2) * sosotb_dim1] = 
2536                                 sosotb[ll + (kk + nd * sosotb_dim2) * 
2537                                 sosotb_dim1] - bid1 * sou1 * sov1 - bid2 * 
2538                                 sou2 * sov1 - bid3 * sou1 * sov2 - bid4 * 
2539                                 sou2 * sov2;
2540                         soditb[ll + (kk + nd * soditb_dim2) * soditb_dim1] = 
2541                                 soditb[ll + (kk + nd * soditb_dim2) * 
2542                                 soditb_dim1] - bid1 * sou1 * div1 - bid2 * 
2543                                 sou2 * div1 - bid3 * sou1 * div2 - bid4 * 
2544                                 sou2 * div2;
2545                         disotb[ll + (kk + nd * disotb_dim2) * disotb_dim1] = 
2546                                 disotb[ll + (kk + nd * disotb_dim2) * 
2547                                 disotb_dim1] - bid1 * diu1 * sov1 - bid2 * 
2548                                 diu2 * sov1 - bid3 * diu1 * sov2 - bid4 * 
2549                                 diu2 * sov2;
2550                         diditb[ll + (kk + nd * diditb_dim2) * diditb_dim1] = 
2551                                 diditb[ll + (kk + nd * diditb_dim2) * 
2552                                 diditb_dim1] - bid1 * diu1 * div1 - bid2 * 
2553                                 diu2 * div1 - bid3 * diu1 * div2 - bid4 * 
2554                                 diu2 * div2;
2555 /* L450: */
2556                     }
2557 /* L400: */
2558                 }
2559
2560 /* ------------ Case when the discretization is done only on the roots  
2561 ----------- */
2562 /* ----------   of Legendre polynom of uneven degree, 0 is root 
2563 ----------- */
2564
2565                 if (*nbpntu % 2 == 1) {
2566                     sou1 = fpntbu[nuroo + 1 + ((ii << 1) - 1) * fpntbu_dim1];
2567                     sou2 = fpntbu[nuroo + 1 + (ii << 1) * fpntbu_dim1];
2568                     i__4 = nvroo;
2569                     for (kk = 1; kk <= i__4; ++kk) {
2570                         kkp = (*nbpntv + 1) / 2 + kk;
2571                         kkm = nvroo - kk + 1;
2572                         sov1 = fpntbv[kkp + ((jj << 1) - 1) * fpntbv_dim1] + 
2573                                 fpntbv[kkm + ((jj << 1) - 1) * fpntbv_dim1];
2574                         div1 = fpntbv[kkp + ((jj << 1) - 1) * fpntbv_dim1] - 
2575                                 fpntbv[kkm + ((jj << 1) - 1) * fpntbv_dim1];
2576                         sov2 = fpntbv[kkp + (jj << 1) * fpntbv_dim1] + fpntbv[
2577                                 kkm + (jj << 1) * fpntbv_dim1];
2578                         div2 = fpntbv[kkp + (jj << 1) * fpntbv_dim1] - fpntbv[
2579                                 kkm + (jj << 1) * fpntbv_dim1];
2580                         sosotb[(kk + nd * sosotb_dim2) * sosotb_dim1] = 
2581                                 sosotb[(kk + nd * sosotb_dim2) * sosotb_dim1] 
2582                                 - bid1 * sou1 * sov1 - bid2 * sou2 * sov1 - 
2583                                 bid3 * sou1 * sov2 - bid4 * sou2 * sov2;
2584                         diditb[(kk + nd * diditb_dim2) * diditb_dim1] = 
2585                                 diditb[(kk + nd * diditb_dim2) * diditb_dim1] 
2586                                 - bid1 * sou1 * div1 - bid2 * sou2 * div1 - 
2587                                 bid3 * sou1 * div2 - bid4 * sou2 * div2;
2588 /* L500: */
2589                     }
2590                 }
2591
2592                 if (*nbpntv % 2 == 1) {
2593                     sov1 = fpntbv[nvroo + 1 + ((jj << 1) - 1) * fpntbv_dim1];
2594                     sov2 = fpntbv[nvroo + 1 + (jj << 1) * fpntbv_dim1];
2595                     i__4 = nuroo;
2596                     for (ll = 1; ll <= i__4; ++ll) {
2597                         llp = (*nbpntu + 1) / 2 + ll;
2598                         llm = nuroo - ll + 1;
2599                         sou1 = fpntbu[llp + ((ii << 1) - 1) * fpntbu_dim1] + 
2600                                 fpntbu[llm + ((ii << 1) - 1) * fpntbu_dim1];
2601                         diu1 = fpntbu[llp + ((ii << 1) - 1) * fpntbu_dim1] - 
2602                                 fpntbu[llm + ((ii << 1) - 1) * fpntbu_dim1];
2603                         sou2 = fpntbu[llp + (ii << 1) * fpntbu_dim1] + fpntbu[
2604                                 llm + (ii << 1) * fpntbu_dim1];
2605                         diu2 = fpntbu[llp + (ii << 1) * fpntbu_dim1] - fpntbu[
2606                                 llm + (ii << 1) * fpntbu_dim1];
2607                         sosotb[ll + nd * sosotb_dim2 * sosotb_dim1] = sosotb[
2608                                 ll + nd * sosotb_dim2 * sosotb_dim1] - bid1 * 
2609                                 sou1 * sov1 - bid2 * sou2 * sov1 - bid3 * 
2610                                 sou1 * sov2 - bid4 * sou2 * sov2;
2611                         diditb[ll + nd * diditb_dim2 * diditb_dim1] = diditb[
2612                                 ll + nd * diditb_dim2 * diditb_dim1] - bid1 * 
2613                                 diu1 * sov1 - bid2 * diu2 * sov1 - bid3 * 
2614                                 diu1 * sov2 - bid4 * diu2 * sov2;
2615 /* L600: */
2616                     }
2617                 }
2618
2619                 if (*nbpntu % 2 == 1 && *nbpntv % 2 == 1) {
2620                     sou1 = fpntbu[nuroo + 1 + ((ii << 1) - 1) * fpntbu_dim1];
2621                     sou2 = fpntbu[nuroo + 1 + (ii << 1) * fpntbu_dim1];
2622                     sov1 = fpntbv[nvroo + 1 + ((jj << 1) - 1) * fpntbv_dim1];
2623                     sov2 = fpntbv[nvroo + 1 + (jj << 1) * fpntbv_dim1];
2624                     sosotb[nd * sosotb_dim2 * sosotb_dim1] = sosotb[nd * 
2625                             sosotb_dim2 * sosotb_dim1] - bid1 * sou1 * sov1 - 
2626                             bid2 * sou2 * sov1 - bid3 * sou1 * sov2 - bid4 * 
2627                             sou2 * sov2;
2628                     diditb[nd * diditb_dim2 * diditb_dim1] = diditb[nd * 
2629                             diditb_dim2 * diditb_dim1] - bid1 * sou1 * sov1 - 
2630                             bid2 * sou2 * sov1 - bid3 * sou1 * sov2 - bid4 * 
2631                             sou2 * sov2;
2632                 }
2633
2634 /* L300: */
2635             }
2636 /* L200: */
2637         }
2638 /* L100: */
2639     }
2640     goto L9999;
2641
2642 /* ------------------------------ The End ------------------------------- 
2643 */
2644
2645 L9999:
2646     if (ibb >= 3) {
2647         AdvApp2Var_SysBase::mgsomsg_("MMA2CD1", 7L);
2648     }
2649     return 0;
2650 } /* mma2cd1_ */
2651
2652 //=======================================================================
2653 //function : mma2cd2_
2654 //purpose  : 
2655 //=======================================================================
2656 int mma2cd2_(integer *ndimen, 
2657              integer *nbpntu, 
2658              integer *nbpntv, 
2659              doublereal *vrootl, 
2660              integer *iordrv, 
2661              doublereal *sotbv1, 
2662              doublereal *sotbv2, 
2663              doublereal *ditbv1, 
2664              doublereal *ditbv2, 
2665              doublereal *fpntab, 
2666              doublereal *vhermt, 
2667              doublereal *sosotb, 
2668              doublereal *soditb, 
2669              doublereal *disotb, 
2670              doublereal *diditb)
2671
2672 {
2673   static integer c__1 = 1;
2674   /* System generated locals */
2675   integer sotbv1_dim1, sotbv1_dim2, sotbv1_offset, sotbv2_dim1, sotbv2_dim2,
2676   sotbv2_offset, ditbv1_dim1, ditbv1_dim2, ditbv1_offset, 
2677   ditbv2_dim1, ditbv2_dim2, ditbv2_offset, fpntab_dim1, 
2678   fpntab_offset, vhermt_dim1, vhermt_offset, sosotb_dim1, 
2679   sosotb_dim2, sosotb_offset, diditb_dim1, diditb_dim2, 
2680   diditb_offset, soditb_dim1, soditb_dim2, soditb_offset, 
2681   disotb_dim1, disotb_dim2, disotb_offset, i__1, i__2, i__3, i__4;
2682
2683   /* Local variables */
2684   static integer ncfhv, nuroo, nvroo, ii, nd, jj, kk, ibb, jjm, jjp;
2685   static doublereal bid1, bid2, bid3, bid4;
2686
2687 /* ********************************************************************** 
2688 */
2689 /*     FUNCTION : */
2690 /*     ---------- */
2691 /*     Discretisation on the parameters of polynoms of interpolation */
2692 /*     of constraints on 2 borders iso-V of order IORDRV. */
2693
2694
2695 /*     KEYWORDS : */
2696 /*     ----------- */
2697 /*     TOUS, AB_SPECIFI::CONTRAINTE&, DISCRETISATION, &POINT */
2698
2699
2700
2701 /*     INPUT ARGUMENTS : */
2702 /*     ------------------ */
2703 /*     NDIMEN: Dimension of the space. */
2704 /*     NBPNTU: Nb of INTERNAL parameters of discretisation by U. */
2705 /*             This is also the nb of root of Legendre polynom where discretization is done. */
2706 /*     UROOTL: Table of parameters of discretisation ON (-1,1) by U. 
2707 */
2708 /*     NBPNTV: Nb of INTERNAL  parameters of discretisation by V. */
2709 /*             This is also the nb of root of Legendre polynom where discretization is done. */
2710 /*     VROOTL: Table of discretization parameters on (-1,1) by V. 
2711 /*     IORDRV: Order of constraint imposed at the extremities of iso-V */
2712 /*             = 0, calculate the extremities of iso-V */
2713 /*             = 1, calculate, additionally, the 1st derivative in the direction of iso-V */
2714 /*             = 2, calculate, additionally, the 2nd derivative in the direction of iso-V */
2715 /*     SOTBV1: Table of NBPNTV/2 sums of 2 index points  */
2716 /*             NBPNTV-II+1 and II, for II = 1, NBPNTV/2 on iso-V0. */
2717 /*     SOTBV2: Table of NBPNTV/2 sums of 2 index points  */
2718 /*             NBPNTV-II+1 and II, for II = 1, NBPNTV/2 on iso-V1. */
2719 /*     DITBV1: Table of NBPNTV/2 differences of 2 index points */
2720 /*             NBPNTV-II+1 and II, for II = 1, NBPNTV/2 on iso-V0. */
2721 /*     DITBV2: Table of NBPNTV/2 differences of 2 index points */
2722 /*             NBPNTV-II+1 and II, for II = 1, NBPNTV/2 on iso-V1. */
2723 /*     SOSOTB: Preinitialized table (input/output argument). */
2724 /*     DISOTB: Preinitialized table (input/output argument). */
2725 /*     SODITB: Preinitialized table (input/output argument). */
2726 /*     DIDITB: Preinitialized table (input/output argument) */
2727
2728 /*     OUTPUT ARGUMENTS : */
2729 /*     ------------------- */
2730 /*     FPNTAB: Auxiliary table. */
2731 /*     VHERMT: Table of 2*(IORDRV+1) coeff. of 2*(IORDRV+1) polynoms of Hermite. */
2732 /*   SOSOTB: Table where the terms of constraints are added */
2733 /*           C(ui,vj) + C(ui,-vj) + C(-ui,vj) + C(-ui,-vj) */
2734 /*           with ui and vj positive roots of the Legendre polynom */
2735 /*           of degree NBPNTU and NBPNTV respectively. */
2736 /*   DISOTB: Table where the terms of constraints are added */
2737 /*           C(ui,vj) + C(ui,-vj) - C(-ui,vj) - C(-ui,-vj) */
2738 /*           with ui and vj positive roots of the polynom of Legendre */
2739 /*           of degree NBPNTU and NBPNTV respectively. */
2740 /*   SODITB: Table where the terms of constraints are added */
2741 /*           C(ui,vj) - C(ui,-vj) + C(-ui,vj) - C(-ui,-vj) */
2742 /*           with ui and vj positive roots of the polynom of Legendre */
2743 /*           of degree NBPNTU and NBPNTV respectively. */
2744 /*   DIDITB: Table where the terms of constraints are added */
2745 /*           C(ui,vj) - C(ui,-vj) - C(-ui,vj) + C(-ui,-vj) */
2746 /*           with ui and vj positive roots of the polynom of Legendre */
2747 /*           of degree NBPNTU and NBPNTV respectively. */
2748
2749 /*     COMMONS USED   : */
2750 /*     ---------------- */
2751
2752 /*     REFERENCES CALLED   : */
2753 /*     ----------------------- */
2754
2755 /*     DESCRIPTION/NOTES/LIMITATIONS : */
2756 /*     ----------------------------------- */
2757
2758
2759 /* > */
2760 /* ********************************************************************** 
2761 */
2762
2763 /*   Name of the routine */
2764
2765
2766     /* Parameter adjustments */
2767     diditb_dim1 = *nbpntu / 2 + 1;
2768     diditb_dim2 = *nbpntv / 2 + 1;
2769     diditb_offset = diditb_dim1 * diditb_dim2;
2770     diditb -= diditb_offset;
2771     disotb_dim1 = *nbpntu / 2;
2772     disotb_dim2 = *nbpntv / 2;
2773     disotb_offset = disotb_dim1 * (disotb_dim2 + 1) + 1;
2774     disotb -= disotb_offset;
2775     soditb_dim1 = *nbpntu / 2;
2776     soditb_dim2 = *nbpntv / 2;
2777     soditb_offset = soditb_dim1 * (soditb_dim2 + 1) + 1;
2778     soditb -= soditb_offset;
2779     sosotb_dim1 = *nbpntu / 2 + 1;
2780     sosotb_dim2 = *nbpntv / 2 + 1;
2781     sosotb_offset = sosotb_dim1 * sosotb_dim2;
2782     sosotb -= sosotb_offset;
2783     --vrootl;
2784     vhermt_dim1 = (*iordrv << 1) + 2;
2785     vhermt_offset = vhermt_dim1;
2786     vhermt -= vhermt_offset;
2787     fpntab_dim1 = *nbpntv;
2788     fpntab_offset = fpntab_dim1 + 1;
2789     fpntab -= fpntab_offset;
2790     ditbv2_dim1 = *nbpntu / 2 + 1;
2791     ditbv2_dim2 = *ndimen;
2792     ditbv2_offset = ditbv2_dim1 * (ditbv2_dim2 + 1);
2793     ditbv2 -= ditbv2_offset;
2794     ditbv1_dim1 = *nbpntu / 2 + 1;
2795     ditbv1_dim2 = *ndimen;
2796     ditbv1_offset = ditbv1_dim1 * (ditbv1_dim2 + 1);
2797     ditbv1 -= ditbv1_offset;
2798     sotbv2_dim1 = *nbpntu / 2 + 1;
2799     sotbv2_dim2 = *ndimen;
2800     sotbv2_offset = sotbv2_dim1 * (sotbv2_dim2 + 1);
2801     sotbv2 -= sotbv2_offset;
2802     sotbv1_dim1 = *nbpntu / 2 + 1;
2803     sotbv1_dim2 = *ndimen;
2804     sotbv1_offset = sotbv1_dim1 * (sotbv1_dim2 + 1);
2805     sotbv1 -= sotbv1_offset;
2806
2807     /* Function Body */
2808     ibb = AdvApp2Var_SysBase::mnfndeb_();
2809     if (ibb >= 3) {
2810         AdvApp2Var_SysBase::mgenmsg_("MMA2CD2", 7L);
2811     }
2812
2813 /* ------------------- Discretization of Hermit polynoms ----------- 
2814 */
2815
2816     ncfhv = (*iordrv + 1) << 1;
2817     i__1 = ncfhv;
2818     for (ii = 1; ii <= i__1; ++ii) {
2819         i__2 = *nbpntv;
2820         for (jj = 1; jj <= i__2; ++jj) {
2821             AdvApp2Var_MathBase::mmmpocur_(&ncfhv, &c__1, &ncfhv, &vhermt[ii * vhermt_dim1], &
2822                     vrootl[jj], &fpntab[jj + ii * fpntab_dim1]);
2823 /* L60: */
2824         }
2825 /* L50: */
2826     }
2827
2828 /* ---- The discretizations of polynoms of constraints are subtracted ---- 
2829 */
2830
2831     nuroo = *nbpntu / 2;
2832     nvroo = *nbpntv / 2;
2833
2834     i__1 = *ndimen;
2835     for (nd = 1; nd <= i__1; ++nd) {
2836         i__2 = *iordrv + 1;
2837         for (ii = 1; ii <= i__2; ++ii) {
2838
2839             i__3 = nuroo;
2840             for (kk = 1; kk <= i__3; ++kk) {
2841                 bid1 = sotbv1[kk + (nd + ii * sotbv1_dim2) * sotbv1_dim1];
2842                 bid2 = sotbv2[kk + (nd + ii * sotbv2_dim2) * sotbv2_dim1];
2843                 bid3 = ditbv1[kk + (nd + ii * ditbv1_dim2) * ditbv1_dim1];
2844                 bid4 = ditbv2[kk + (nd + ii * ditbv2_dim2) * ditbv2_dim1];
2845                 i__4 = nvroo;
2846                 for (jj = 1; jj <= i__4; ++jj) {
2847                     jjp = (*nbpntv + 1) / 2 + jj;
2848                     jjm = nvroo - jj + 1;
2849                     sosotb[kk + (jj + nd * sosotb_dim2) * sosotb_dim1] = 
2850                             sosotb[kk + (jj + nd * sosotb_dim2) * sosotb_dim1]
2851                              - bid1 * (fpntab[jjp + ((ii << 1) - 1) * 
2852                             fpntab_dim1] + fpntab[jjm + ((ii << 1) - 1) * 
2853                             fpntab_dim1]) - bid2 * (fpntab[jjp + (ii << 1) * 
2854                             fpntab_dim1] + fpntab[jjm + (ii << 1) * 
2855                             fpntab_dim1]);
2856                     disotb[kk + (jj + nd * disotb_dim2) * disotb_dim1] = 
2857                             disotb[kk + (jj + nd * disotb_dim2) * disotb_dim1]
2858                              - bid3 * (fpntab[jjp + ((ii << 1) - 1) * 
2859                             fpntab_dim1] + fpntab[jjm + ((ii << 1) - 1) * 
2860                             fpntab_dim1]) - bid4 * (fpntab[jjp + (ii << 1) * 
2861                             fpntab_dim1] + fpntab[jjm + (ii << 1) * 
2862                             fpntab_dim1]);
2863                     soditb[kk + (jj + nd * soditb_dim2) * soditb_dim1] = 
2864                             soditb[kk + (jj + nd * soditb_dim2) * soditb_dim1]
2865                              - bid1 * (fpntab[jjp + ((ii << 1) - 1) * 
2866                             fpntab_dim1] - fpntab[jjm + ((ii << 1) - 1) * 
2867                             fpntab_dim1]) - bid2 * (fpntab[jjp + (ii << 1) * 
2868                             fpntab_dim1] - fpntab[jjm + (ii << 1) * 
2869                             fpntab_dim1]);
2870                     diditb[kk + (jj + nd * diditb_dim2) * diditb_dim1] = 
2871                             diditb[kk + (jj + nd * diditb_dim2) * diditb_dim1]
2872                              - bid3 * (fpntab[jjp + ((ii << 1) - 1) * 
2873                             fpntab_dim1] - fpntab[jjm + ((ii << 1) - 1) * 
2874                             fpntab_dim1]) - bid4 * (fpntab[jjp + (ii << 1) * 
2875                             fpntab_dim1] - fpntab[jjm + (ii << 1) * 
2876                             fpntab_dim1]);
2877 /* L400: */
2878                 }
2879 /* L300: */
2880             }
2881 /* L200: */
2882         }
2883
2884 /* ------------ Case when the discretization is done only on the roots  */
2885 /* ----------   of Legendre polynom of uneven degree, 0 is root */
2886
2887
2888         if (*nbpntv % 2 == 1) {
2889             i__2 = *iordrv + 1;
2890             for (ii = 1; ii <= i__2; ++ii) {
2891                 i__3 = nuroo;
2892                 for (kk = 1; kk <= i__3; ++kk) {
2893                     bid1 = sotbv1[kk + (nd + ii * sotbv1_dim2) * sotbv1_dim1] 
2894                             * fpntab[nvroo + 1 + ((ii << 1) - 1) * 
2895                             fpntab_dim1] + sotbv2[kk + (nd + ii * sotbv2_dim2)
2896                              * sotbv2_dim1] * fpntab[nvroo + 1 + (ii << 1) * 
2897                             fpntab_dim1];
2898                     sosotb[kk + nd * sosotb_dim2 * sosotb_dim1] -= bid1;
2899                     bid2 = ditbv1[kk + (nd + ii * ditbv1_dim2) * ditbv1_dim1] 
2900                             * fpntab[nvroo + 1 + ((ii << 1) - 1) * 
2901                             fpntab_dim1] + ditbv2[kk + (nd + ii * ditbv2_dim2)
2902                              * ditbv2_dim1] * fpntab[nvroo + 1 + (ii << 1) * 
2903                             fpntab_dim1];
2904                     diditb[kk + nd * diditb_dim2 * diditb_dim1] -= bid2;
2905 /* L550: */
2906                 }
2907 /* L500: */
2908             }
2909         }
2910
2911         if (*nbpntu % 2 == 1) {
2912             i__2 = *iordrv + 1;
2913             for (ii = 1; ii <= i__2; ++ii) {
2914                 i__3 = nvroo;
2915                 for (jj = 1; jj <= i__3; ++jj) {
2916                     jjp = (*nbpntv + 1) / 2 + jj;
2917                     jjm = nvroo - jj + 1;
2918                     bid1 = sotbv1[(nd + ii * sotbv1_dim2) * sotbv1_dim1] * (
2919                             fpntab[jjp + ((ii << 1) - 1) * fpntab_dim1] + 
2920                             fpntab[jjm + ((ii << 1) - 1) * fpntab_dim1]) + 
2921                             sotbv2[(nd + ii * sotbv2_dim2) * sotbv2_dim1] * (
2922                             fpntab[jjp + (ii << 1) * fpntab_dim1] + fpntab[
2923                             jjm + (ii << 1) * fpntab_dim1]);
2924                     sosotb[(jj + nd * sosotb_dim2) * sosotb_dim1] -= bid1;
2925                     bid2 = sotbv1[(nd + ii * sotbv1_dim2) * sotbv1_dim1] * (
2926                             fpntab[jjp + ((ii << 1) - 1) * fpntab_dim1] - 
2927                             fpntab[jjm + ((ii << 1) - 1) * fpntab_dim1]) + 
2928                             sotbv2[(nd + ii * sotbv2_dim2) * sotbv2_dim1] * (
2929                             fpntab[jjp + (ii << 1) * fpntab_dim1] - fpntab[
2930                             jjm + (ii << 1) * fpntab_dim1]);
2931                     diditb[jj + nd * diditb_dim2 * diditb_dim1] -= bid2;
2932 /* L650: */
2933                 }
2934 /* L600: */
2935             }
2936         }
2937
2938         if (*nbpntu % 2 == 1 && *nbpntv % 2 == 1) {
2939             i__2 = *iordrv + 1;
2940             for (ii = 1; ii <= i__2; ++ii) {
2941                 bid1 = sotbv1[(nd + ii * sotbv1_dim2) * sotbv1_dim1] * fpntab[
2942                         nvroo + 1 + ((ii << 1) - 1) * fpntab_dim1] + sotbv2[(
2943                         nd + ii * sotbv2_dim2) * sotbv2_dim1] * fpntab[nvroo 
2944                         + 1 + (ii << 1) * fpntab_dim1];
2945                 sosotb[nd * sosotb_dim2 * sosotb_dim1] -= bid1;
2946 /* L700: */
2947             }
2948         }
2949
2950 /* L100: */
2951     }
2952     goto L9999;
2953
2954 /* ------------------------------ The End ------------------------------- 
2955 */
2956
2957 L9999:
2958     if (ibb >= 3) {
2959         AdvApp2Var_SysBase::mgsomsg_("MMA2CD2", 7L);
2960     }
2961     return 0;
2962 } /* mma2cd2_ */
2963
2964 //=======================================================================
2965 //function : mma2cd3_
2966 //purpose  : 
2967 //=======================================================================
2968 int mma2cd3_(integer *ndimen,
2969              integer *nbpntu, 
2970              doublereal *urootl, 
2971              integer *nbpntv, 
2972              integer *iordru, 
2973              doublereal *sotbu1, 
2974              doublereal *sotbu2, 
2975              doublereal *ditbu1, 
2976              doublereal *ditbu2, 
2977              doublereal *fpntab, 
2978              doublereal *uhermt, 
2979              doublereal *sosotb, 
2980              doublereal *soditb, 
2981              doublereal *disotb, 
2982              doublereal *diditb)
2983
2984 {
2985   static integer c__1 = 1;
2986
2987    /* System generated locals */
2988     integer sotbu1_dim1, sotbu1_dim2, sotbu1_offset, sotbu2_dim1, sotbu2_dim2,
2989              sotbu2_offset, ditbu1_dim1, ditbu1_dim2, ditbu1_offset, 
2990             ditbu2_dim1, ditbu2_dim2, ditbu2_offset, fpntab_dim1, 
2991             fpntab_offset, uhermt_dim1, uhermt_offset, sosotb_dim1, 
2992             sosotb_dim2, sosotb_offset, diditb_dim1, diditb_dim2, 
2993             diditb_offset, soditb_dim1, soditb_dim2, soditb_offset, 
2994             disotb_dim1, disotb_dim2, disotb_offset, i__1, i__2, i__3, i__4;
2995
2996     /* Local variables */
2997     static integer ncfhu, nuroo, nvroo, ii, nd, jj, kk, ibb, kkm, kkp;
2998     static doublereal bid1, bid2, bid3, bid4;
2999
3000 /* ********************************************************************** 
3001 */
3002 /*     FUNCTION : */
3003 /*     ---------- */
3004 /*     Discretisation on the parameters of polynoms of interpolation */
3005 /*     of constraints on 2 borders iso-U of order IORDRU. */
3006
3007
3008 /*     KEYWORDS : */
3009 /*     ----------- */
3010 /*     TOUS, AB_SPECIFI::CONTRAINTE&, DISCRETISATION, &POINT */
3011
3012 /*     INPUT ARGUMENTS : */
3013 /*     ------------------ */
3014 /*     NDIMEN: Dimension of the space. */
3015 /*     NBPNTU: Nb of INTERNAL parameters of discretisation by U. */
3016 /*             This is also the nb of root of Legendre polynom where discretization is done. */
3017 /*     UROOTL: Table of parameters of discretisation ON (-1,1) by U. 
3018 */
3019 /*     NBPNTV: Nb of INTERNAL  parameters of discretisation by V. */
3020 /*             This is also the nb of root of Legendre polynom where discretization is done. */
3021 /*     IORDRV: Order of constraint imposed at the extremities of iso-V */
3022 /*             = 0, calculate the extremities of iso-V */
3023 /*             = 1, calculate, additionally, the 1st derivative in the direction of iso-V */
3024 /*             = 2, calculate, additionally, the 2nd derivative in the direction of iso-V */
3025 /*     SOTBU1: Table of NBPNTU/2 sums of 2 index points  */
3026 /*             NBPNTU-II+1 and II, for II = 1, NBPNTU/2 on iso-V0. */
3027 /*     SOTBU2: Table of NBPNTV/2 sums of 2 index points  */
3028 /*             NBPNTU-II+1 and II, for II = 1, NBPNTU/2 on iso-V1. */
3029 /*     DITBU1: Table of NBPNTU/2 differences of 2 index points */
3030 /*             NBPNTU-II+1 and II, for II = 1, NBPNTU/2 on iso-V0. */
3031 /*     DITBU2: Table of NBPNTU/2 differences of 2 index points */
3032 /*             NBPNTU-II+1 and II, for II = 1, NBPNTU/2 on iso-V1. */
3033 /*     SOSOTB: Preinitialized table (input/output argument). */
3034 /*     DISOTB: Preinitialized table (input/output argument). */
3035 /*     SODITB: Preinitialized table (input/output argument). */
3036 /*     DIDITB: Preinitialized table (input/output argument) */
3037
3038 /*     OUTPUT ARGUMENTS : */
3039 /*     ------------------- */
3040 /*     FPNTAB: Auxiliary table. */
3041 /*     UHERMT: Table of 2*(IORDRU+1) coeff. of 2*(IORDRU+1) polynoms of Hermite. */
3042 /*   SOSOTB: Table where the terms of constraints are added */
3043 /*           C(ui,vj) + C(ui,-vj) + C(-ui,vj) + C(-ui,-vj) */
3044 /*           with ui and vj positive roots of the Legendre polynom */
3045 /*           of degree NBPNTU and NBPNTV respectively. */
3046 /*   DISOTB: Table where the terms of constraints are added */
3047 /*           C(ui,vj) + C(ui,-vj) - C(-ui,vj) - C(-ui,-vj) */
3048 /*           with ui and vj positive roots of the polynom of Legendre */
3049 /*           of degree NBPNTU and NBPNTV respectively. */
3050 /*   SODITB: Table where the terms of constraints are added */
3051 /*           C(ui,vj) - C(ui,-vj) + C(-ui,vj) - C(-ui,-vj) */
3052 /*           with ui and vj positive roots of the polynom of Legendre */
3053 /*           of degree NBPNTU and NBPNTV respectively. */
3054 /*   DIDITB: Table where the terms of constraints are added */
3055 /*           C(ui,vj) - C(ui,-vj) - C(-ui,vj) + C(-ui,-vj) */
3056 /*           with ui and vj positive roots of the polynom of Legendre */
3057 /*           of degree NBPNTU and NBPNTV respectively. */
3058
3059 /*     COMMONS USED   : */
3060 /*     ---------------- */
3061
3062 /*     REFERENCES CALLED   : */
3063 /*     ----------------------- */
3064
3065 /*     DESCRIPTION/NOTES/LIMITATIONS : */
3066 /*     ----------------------------------- */
3067
3068 /* $    HISTORIQUE DES MODIFICATIONS   : */
3069 /*     -------------------------------- */
3070 /*     08-08-1991: RBD; Creation. */
3071 /* > */
3072 /* ********************************************************************** 
3073 */
3074
3075 /*   Name of the routine */
3076
3077
3078     /* Parameter adjustments */
3079     --urootl;
3080     diditb_dim1 = *nbpntu / 2 + 1;
3081     diditb_dim2 = *nbpntv / 2 + 1;
3082     diditb_offset = diditb_dim1 * diditb_dim2;
3083     diditb -= diditb_offset;
3084     disotb_dim1 = *nbpntu / 2;
3085     disotb_dim2 = *nbpntv / 2;
3086     disotb_offset = disotb_dim1 * (disotb_dim2 + 1) + 1;
3087     disotb -= disotb_offset;
3088     soditb_dim1 = *nbpntu / 2;
3089     soditb_dim2 = *nbpntv / 2;
3090     soditb_offset = soditb_dim1 * (soditb_dim2 + 1) + 1;
3091     soditb -= soditb_offset;
3092     sosotb_dim1 = *nbpntu / 2 + 1;
3093     sosotb_dim2 = *nbpntv / 2 + 1;
3094     sosotb_offset = sosotb_dim1 * sosotb_dim2;
3095     sosotb -= sosotb_offset;
3096     uhermt_dim1 = (*iordru << 1) + 2;
3097     uhermt_offset = uhermt_dim1;
3098     uhermt -= uhermt_offset;
3099     fpntab_dim1 = *nbpntu;
3100     fpntab_offset = fpntab_dim1 + 1;
3101     fpntab -= fpntab_offset;
3102     ditbu2_dim1 = *nbpntv / 2 + 1;
3103     ditbu2_dim2 = *ndimen;
3104     ditbu2_offset = ditbu2_dim1 * (ditbu2_dim2 + 1);
3105     ditbu2 -= ditbu2_offset;
3106     ditbu1_dim1 = *nbpntv / 2 + 1;
3107     ditbu1_dim2 = *ndimen;
3108     ditbu1_offset = ditbu1_dim1 * (ditbu1_dim2 + 1);
3109     ditbu1 -= ditbu1_offset;
3110     sotbu2_dim1 = *nbpntv / 2 + 1;
3111     sotbu2_dim2 = *ndimen;
3112     sotbu2_offset = sotbu2_dim1 * (sotbu2_dim2 + 1);
3113     sotbu2 -= sotbu2_offset;
3114     sotbu1_dim1 = *nbpntv / 2 + 1;
3115     sotbu1_dim2 = *ndimen;
3116     sotbu1_offset = sotbu1_dim1 * (sotbu1_dim2 + 1);
3117     sotbu1 -= sotbu1_offset;
3118
3119     /* Function Body */
3120     ibb = AdvApp2Var_SysBase::mnfndeb_();
3121     if (ibb >= 3) {
3122         AdvApp2Var_SysBase::mgenmsg_("MMA2CD3", 7L);
3123     }
3124
3125 /* ------------------- Discretization of polynoms of Hermit ----------- 
3126 */
3127
3128     ncfhu = (*iordru + 1) << 1;
3129     i__1 = ncfhu;
3130     for (ii = 1; ii <= i__1; ++ii) {
3131         i__2 = *nbpntu;
3132         for (kk = 1; kk <= i__2; ++kk) {
3133             AdvApp2Var_MathBase::mmmpocur_(&ncfhu, 
3134                                            &c__1, 
3135                                            &ncfhu,
3136                                            &uhermt[ii * uhermt_dim1],
3137                                            &urootl[kk], 
3138                                            &fpntab[kk + ii * fpntab_dim1]);
3139 /* L60: */
3140         }
3141 /* L50: */
3142     }
3143
3144 /* ---- The discretizations of polynoms of constraints are subtracted ---- 
3145 */
3146
3147     nvroo = *nbpntv / 2;
3148     nuroo = *nbpntu / 2;
3149
3150     i__1 = *ndimen;
3151     for (nd = 1; nd <= i__1; ++nd) {
3152         i__2 = *iordru + 1;
3153         for (ii = 1; ii <= i__2; ++ii) {
3154
3155             i__3 = nvroo;
3156             for (jj = 1; jj <= i__3; ++jj) {
3157                 bid1 = sotbu1[jj + (nd + ii * sotbu1_dim2) * sotbu1_dim1];
3158                 bid2 = sotbu2[jj + (nd + ii * sotbu2_dim2) * sotbu2_dim1];
3159                 bid3 = ditbu1[jj + (nd + ii * ditbu1_dim2) * ditbu1_dim1];
3160                 bid4 = ditbu2[jj + (nd + ii * ditbu2_dim2) * ditbu2_dim1];
3161                 i__4 = nuroo;
3162                 for (kk = 1; kk <= i__4; ++kk) {
3163                     kkp = (*nbpntu + 1) / 2 + kk;
3164                     kkm = nuroo - kk + 1;
3165                     sosotb[kk + (jj + nd * sosotb_dim2) * sosotb_dim1] = 
3166                             sosotb[kk + (jj + nd * sosotb_dim2) * sosotb_dim1]
3167                              - bid1 * (fpntab[kkp + ((ii << 1) - 1) * 
3168                             fpntab_dim1] + fpntab[kkm + ((ii << 1) - 1) * 
3169                             fpntab_dim1]) - bid2 * (fpntab[kkp + (ii << 1) * 
3170                             fpntab_dim1] + fpntab[kkm + (ii << 1) * 
3171                             fpntab_dim1]);
3172                     disotb[kk + (jj + nd * disotb_dim2) * disotb_dim1] = 
3173                             disotb[kk + (jj + nd * disotb_dim2) * disotb_dim1]
3174                              - bid1 * (fpntab[kkp + ((ii << 1) - 1) * 
3175                             fpntab_dim1] - fpntab[kkm + ((ii << 1) - 1) * 
3176                             fpntab_dim1]) - bid2 * (fpntab[kkp + (ii << 1) * 
3177                             fpntab_dim1] - fpntab[kkm + (ii << 1) * 
3178                             fpntab_dim1]);
3179                     soditb[kk + (jj + nd * soditb_dim2) * soditb_dim1] = 
3180                             soditb[kk + (jj + nd * soditb_dim2) * soditb_dim1]
3181                              - bid3 * (fpntab[kkp + ((ii << 1) - 1) * 
3182                             fpntab_dim1] + fpntab[kkm + ((ii << 1) - 1) * 
3183                             fpntab_dim1]) - bid4 * (fpntab[kkp + (ii << 1) * 
3184                             fpntab_dim1] + fpntab[kkm + (ii << 1) * 
3185                             fpntab_dim1]);
3186                     diditb[kk + (jj + nd * diditb_dim2) * diditb_dim1] = 
3187                             diditb[kk + (jj + nd * diditb_dim2) * diditb_dim1]
3188                              - bid3 * (fpntab[kkp + ((ii << 1) - 1) * 
3189                             fpntab_dim1] - fpntab[kkm + ((ii << 1) - 1) * 
3190                             fpntab_dim1]) - bid4 * (fpntab[kkp + (ii << 1) * 
3191                             fpntab_dim1] - fpntab[kkm + (ii << 1) * 
3192                             fpntab_dim1]);
3193 /* L400: */
3194                 }
3195 /* L300: */
3196             }
3197 /* L200: */
3198         }
3199
3200 /* ------------ Case when the discretization is done only on the roots  */
3201 /* ----------   of Legendre polynom of uneven degree, 0 is root */
3202
3203
3204
3205         if (*nbpntu % 2 == 1) {
3206             i__2 = *iordru + 1;
3207             for (ii = 1; ii <= i__2; ++ii) {
3208                 i__3 = nvroo;
3209                 for (jj = 1; jj <= i__3; ++jj) {
3210                     bid1 = sotbu1[jj + (nd + ii * sotbu1_dim2) * sotbu1_dim1] 
3211                             * fpntab[nuroo + 1 + ((ii << 1) - 1) * 
3212                             fpntab_dim1] + sotbu2[jj + (nd + ii * sotbu2_dim2)
3213                              * sotbu2_dim1] * fpntab[nuroo + 1 + (ii << 1) * 
3214                             fpntab_dim1];
3215                     sosotb[(jj + nd * sosotb_dim2) * sosotb_dim1] -= bid1;
3216 &nbs