fc1cb110c6cf0d760e2576b56d6121b2006772c1
[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     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                 foncnp.Evaluate (ndimen, &uvfonc[3], &uvfonc[5], isofav, tconst,
874                         &nbp, ttable, &ideru, &iderv, &contr1[(iderv + 1) * 
875                         contr1_dim1 + 1], iercod);
876                 if (*iercod > 0) {
877                     goto L9999;
878                 }
879 /* L500: */
880             }
881             i__1 = *iordre;
882             for (iderv = 1; iderv <= i__1; ++iderv) {
883                 foncnp.Evaluate (ndimen, &uvfonc[3], &uvfonc[5], isofav, tconst,
884                         &nbp, &ttable[*nbroot + 1], &ideru, &iderv, &contr2[(
885                         iderv + 1) * contr2_dim1 + 1], iercod);
886                 if (*iercod > 0) {
887                     goto L9999;
888                 }
889 /* L510: */
890             }
891 /*    If Iso-V, derive by U till order IORDRE */
892         } else {
893 /* --> Factor of normalization  1st derivative. */
894             bid1 = (uvfonc[4] - uvfonc[3]) / 2.;
895             i__1 = *iordre;
896             for (ideru = 1; ideru <= i__1; ++ideru) {
897                 foncnp.Evaluate (ndimen, &uvfonc[3], &uvfonc[5], isofav, tconst,
898                         &nbp, ttable, &ideru, &iderv, &contr1[(ideru + 1) * 
899                         contr1_dim1 + 1], iercod);
900                 if (*iercod > 0) {
901                     goto L9999;
902                 }
903 /* L600: */
904             }
905             i__1 = *iordre;
906             for (ideru = 1; ideru <= i__1; ++ideru) {
907                 foncnp.Evaluate (ndimen, &uvfonc[3], &uvfonc[5], isofav, tconst,
908                         &nbp, &ttable[*nbroot + 1], &ideru, &iderv, &contr2[(
909                         ideru + 1) * contr2_dim1 + 1], iercod);
910                 if (*iercod > 0) {
911                     goto L9999;
912                 }
913 /* L610: */
914             }
915         }
916
917 /* ------------------------- Normalization of derivatives -------------
918 ---- */
919 /* (The function is redefined on (-1,1)*(-1,1)) */
920         bid2 = renor;
921         i__1 = *iordre;
922         for (ii = 1; ii <= i__1; ++ii) {
923             bid2 = bid1 * bid2;
924             i__2 = *ndimen;
925             for (nd = 1; nd <= i__2; ++nd) {
926                 contr1[nd + (ii + 1) * contr1_dim1] *= bid2;
927                 contr2[nd + (ii + 1) * contr2_dim1] *= bid2;
928 /* L710: */
929             }
930 /* L700: */
931         }
932     }
933
934 /* ------------------------------ The end ------------------------------- 
935 */
936
937 L9999:
938     if (*iercod > 0) {
939         *iercod += 100;
940         AdvApp2Var_SysBase::maermsg_("MMA1FDI", iercod, 7L);
941     }
942     if (ibb >= 3) {
943         AdvApp2Var_SysBase::mgsomsg_("MMA1FDI", 7L);
944     }
945     return 0;
946 } /* mma1fdi_ */
947
948 //=======================================================================
949 //function : mma1fer_
950 //purpose  : 
951 //=======================================================================
952 int mma1fer_(integer *,//ndimen, 
953              integer *nbsesp, 
954              integer *ndimse, 
955              integer *iordre, 
956              integer *ndgjac, 
957              doublereal *crvjac, 
958              integer *ncflim, 
959              doublereal *epsapr, 
960              doublereal *ycvmax, 
961              doublereal *errmax, 
962              doublereal *errmoy, 
963              integer *ncoeff, 
964              integer *iercod)
965 {
966   /* System generated locals */
967   integer crvjac_dim1, crvjac_offset, i__1, i__2;
968
969   /* Local variables */
970   static integer idim, ncfja, ncfnw, ndses, ii, kk, ibb, ier;
971   static integer nbr0;
972
973
974 /* ***********************************************************************
975  */
976
977 /*     FUNCTION : */
978 /*     ---------- */
979 /*     Calculate the degree and the errors of approximation of a border. */
980
981 /*     KEYWORDS : */
982 /*     ----------- */
983 /*      TOUS,AB_SPECIFI :: COURBE&,TRONCATURE, &PRECISION */
984
985 /*     INPUT ARGUMENTS : */
986 /*     -------------------- */
987
988 /*     NDIMEN: Total Dimension of the space (sum of dimensions of sub-spaces) */
989 /*     NBSESP: Number of "independent" sub-spaces. */
990 /*     NDIMSE: Table of dimensions of sub-spaces. */
991 /*     IORDRE: Order of constraint at the extremities of the border */
992 /*              -1 = no constraints, */
993 /*               0 = constraints of passage to limits (i.e. C0), */
994 /*               1 = C0 + constraintes of 1st derivatives (i.e. C1), */
995 /*               2 = C1 + constraintes of 2nd derivatives (i.e. C2). */
996 /*     NDGJAC: Degree of development in series to use for the calculation  
997 /*             in the base of Jacobi. */
998 /*     CRVJAC: Table of coeff. of the curve of approximation in the */
999 /*             base of Jacobi. */
1000 /*     NCFLIM: Max number of coeff of the polynomial curve */
1001 /*             of approximation (should be above or equal to */
1002 /*             2*IORDRE+2 and below or equal to 50). */
1003 /*     EPSAPR: Table of errors of approximations that cannot be passed, */
1004 /*             sub-space by sub-space. */
1005
1006 /*     OUTPUT ARGUMENTS : */
1007 /*     --------------------- */
1008 /*     YCVMAX: Auxiliary Table. */
1009 /*     ERRMAX: Table of errors (sub-space by sub-space) */
1010 /*             MAXIMUM made in the approximation of FONCNP by */
1011 /*             COURBE. */
1012 /*     ERRMOY: Table of errors (sub-space by sub-space) */
1013 /*             AVERAGE made in the approximation of FONCNP by */
1014 /*             COURBE. */
1015 /*     NCOEFF: Number of significative coeffs. of the calculated "curve". */
1016 /*     IERCOD: Error code */
1017 /*             = 0, ok, */
1018 /*             =-1, warning, required tolerance can't be */
1019 /*                  met with coefficients NFCLIM. */
1020 /*             = 1, order of constraints (IORDRE) is not within authorised values */
1021 /*                  
1022
1023 /*     COMMONS USED : */
1024 /*     ------------------ */
1025
1026 /*     REFERENCES CALLED : */
1027 /*     --------------------- */
1028
1029 /*     DESCRIPTION/NOTES/LIMITATIONS : */
1030 /*     ----------------------------------- */
1031 /* > */
1032 /* ********************************************************************** 
1033 */
1034
1035 /*  Name of the routine */
1036
1037
1038     /* Parameter adjustments */
1039     --ycvmax;
1040     --errmoy;
1041     --errmax;
1042     --epsapr;
1043     --ndimse;
1044     crvjac_dim1 = *ndgjac + 1;
1045     crvjac_offset = crvjac_dim1;
1046     crvjac -= crvjac_offset;
1047
1048     /* Function Body */
1049     ibb = AdvApp2Var_SysBase::mnfndeb_();
1050     if (ibb >= 3) {
1051         AdvApp2Var_SysBase::mgenmsg_("MMA1FER", 7L);
1052     }
1053     *iercod = 0;
1054     idim = 1;
1055     *ncoeff = 0;
1056     ncfja = *ndgjac + 1;
1057
1058 /* ------------ Calculate the degree of the curve and of the Max error -------- 
1059 */
1060 /* -------------- of approximation for all sub-spaces -------- 
1061 */
1062
1063     i__1 = *nbsesp;
1064     for (ii = 1; ii <= i__1; ++ii) {
1065         ndses = ndimse[ii];
1066
1067 /* ------------ cutting of coeff. and calculation of Max error -------
1068 ---- */
1069
1070         AdvApp2Var_MathBase::mmtrpjj_(&ncfja, &ndses, &ncfja, &epsapr[ii], iordre, &crvjac[idim * 
1071                 crvjac_dim1], &ycvmax[1], &errmax[ii], &ncfnw);
1072
1073 /* ******************************************************************
1074 **** */
1075 /* ------------- If precision OK, calculate the average error -------
1076 ---- */
1077 /* ******************************************************************
1078 **** */
1079
1080         if (ncfnw <= *ncflim) {
1081             mmaperm_(&ncfja, &ndses, &ncfja, iordre, &crvjac[idim * 
1082                     crvjac_dim1], &ncfnw, &errmoy[ii]);
1083             *ncoeff = advapp_max(ncfnw,*ncoeff);
1084
1085 /* ------------- Set the declined coefficients to 0.D0 -----------
1086 -------- */
1087
1088             nbr0 = *ncflim - ncfnw;
1089             if (nbr0 > 0) {
1090                 i__2 = ndses;
1091                 for (kk = 1; kk <= i__2; ++kk) {
1092                   AdvApp2Var_SysBase::mvriraz_(&nbr0, 
1093                              (char *)&crvjac[ncfnw + (idim + kk - 1) * crvjac_dim1]);
1094 /* L200: */
1095                 }
1096             }
1097         } else {
1098
1099 /* **************************************************************
1100 ******** */
1101 /* ------------------- If required precision can't be reached----
1102 -------- */
1103 /* **************************************************************
1104 ******** */
1105
1106             *iercod = -1;
1107
1108 /* ------------------------- calculate the Max error ------------
1109 -------- */
1110
1111             AdvApp2Var_MathBase::mmaperx_(&ncfja, &ndses, &ncfja, iordre, &crvjac[idim * 
1112                     crvjac_dim1], ncflim, &ycvmax[1], &errmax[ii], &ier);
1113             if (ier > 0) {
1114                 goto L9100;
1115             }
1116
1117 /* -------------------- nb of coeff to be returned -------------
1118 -------- */
1119
1120             *ncoeff = *ncflim;
1121
1122 /* ------------------- and calculation of the average error ----
1123 -------- */
1124
1125             mmaperm_(&ncfja, &ndses, &ncfja, iordre, &crvjac[idim * 
1126                     crvjac_dim1], ncflim, &errmoy[ii]);
1127         }
1128         idim += ndses;
1129 /* L100: */
1130     }
1131
1132     goto L9999;
1133
1134 /* ------------------------------ The end ------------------------------- 
1135 */
1136 /* --> The order of constraints is not within autorized values. */
1137 L9100:
1138     *iercod = 1;
1139     goto L9999;
1140
1141 L9999:
1142     if (*iercod != 0) {
1143         AdvApp2Var_SysBase::maermsg_("MMA1FER", iercod, 7L);
1144     }
1145     if (ibb >= 3) {
1146         AdvApp2Var_SysBase::mgsomsg_("MMA1FER", 7L);
1147     }
1148     return 0;
1149 } /* mma1fer_ */
1150
1151
1152 //=======================================================================
1153 //function : mma1her_
1154 //purpose  : 
1155 //=======================================================================
1156 int AdvApp2Var_ApproxF2var::mma1her_(const integer *iordre, 
1157                                      doublereal *hermit, 
1158                                      integer *iercod)
1159 {
1160   /* System generated locals */
1161   integer hermit_dim1, hermit_offset;
1162
1163   /* Local variables */
1164   static integer ibb;
1165
1166
1167
1168 /* ********************************************************************** 
1169 */
1170
1171 /*     FUNCTION : */
1172 /*     ---------- */
1173 /*     Calculate 2*(IORDRE+1) Hermit polynoms of  degree 2*IORDRE+1 */
1174 /*     on (-1,1) */
1175
1176 /*     KEYWORDS : */
1177 /*     ----------- */
1178 /*     ALL, AB_SPECIFI::CONTRAINTE&, INTERPOLATION, &POLYNOME */
1179
1180 /*     INPUT ARGUMENTS : */
1181 /*     ------------------ */
1182 /*     IORDRE: Order of constraint. */
1183 /*      = 0, Polynom of interpolation of order C0 on (-1,1). */
1184 /*      = 1, Polynom of interpolation of order C0 and C1 on (-1,1). */
1185 /*      = 2, Polynom of interpolation of order C0, C1 and C2 on (-1,1). 
1186 */
1187
1188 /*     OUTPUT ARGUMENTS : */
1189 /*     ------------------- */
1190 /*     HERMIT: Table of 2*IORDRE+2 coeff. of each of  2*(IORDRE+1) */
1191 /*             HERMIT polynom. */
1192 /*     IERCOD: Error code, */
1193 /*      = 0, Ok */
1194 /*      = 1, required order of constraint is not managed here. */
1195 /*     COMMONS USED   : */
1196 /*     ---------------- */
1197
1198 /*     REFERENCES CALLED   : */
1199 /*     ----------------------- */
1200
1201 /*     DESCRIPTION/NOTES/LIMITATIONS : */
1202 /*     ----------------------------------- */
1203 /*     The part of HERMIT(*,2*i+j) table where  j=1 or 2 and i=0 to IORDRE, 
1204 /*     contains the coefficients of the polynom of degree 2*IORDRE+1 */
1205 /*     such as ALL values in -1 and in +1 of this polynom and its */
1206 /*     derivatives till order of derivation IORDRE are NULL, */
1207 /*     EXCEPT for the derivative of order i: */
1208 /*        - valued 1 in -1 if j=1 */
1209 /*        - valued 1 in +1 if j=2. */
1210 /* > */
1211 /* ********************************************************************** 
1212 */
1213
1214 /*  Name of the routine */
1215
1216
1217     /* Parameter adjustments */
1218     hermit_dim1 = (*iordre + 1) << 1;
1219     hermit_offset = hermit_dim1 + 1;
1220     hermit -= hermit_offset;
1221
1222     /* Function Body */
1223     ibb = AdvApp2Var_SysBase::mnfndeb_();
1224     if (ibb >= 3) {
1225         AdvApp2Var_SysBase::mgenmsg_("MMA1HER", 7L);
1226     }
1227     *iercod = 0;
1228
1229 /* --- Recover (IORDRE+2) coeff of 2*(IORDRE+1) Hermit polynoms -- 
1230 */
1231
1232     if (*iordre == 0) {
1233         hermit[hermit_dim1 + 1] = .5;
1234         hermit[hermit_dim1 + 2] = -.5;
1235
1236         hermit[(hermit_dim1 << 1) + 1] = .5;
1237         hermit[(hermit_dim1 << 1) + 2] = .5;
1238     } else if (*iordre == 1) {
1239         hermit[hermit_dim1 + 1] = .5;
1240         hermit[hermit_dim1 + 2] = -.75;
1241         hermit[hermit_dim1 + 3] = 0.;
1242         hermit[hermit_dim1 + 4] = .25;
1243
1244         hermit[(hermit_dim1 << 1) + 1] = .5;
1245         hermit[(hermit_dim1 << 1) + 2] = .75;
1246         hermit[(hermit_dim1 << 1) + 3] = 0.;
1247         hermit[(hermit_dim1 << 1) + 4] = -.25;
1248
1249         hermit[hermit_dim1 * 3 + 1] = .25;
1250         hermit[hermit_dim1 * 3 + 2] = -.25;
1251         hermit[hermit_dim1 * 3 + 3] = -.25;
1252         hermit[hermit_dim1 * 3 + 4] = .25;
1253
1254         hermit[(hermit_dim1 << 2) + 1] = -.25;
1255         hermit[(hermit_dim1 << 2) + 2] = -.25;
1256         hermit[(hermit_dim1 << 2) + 3] = .25;
1257         hermit[(hermit_dim1 << 2) + 4] = .25;
1258     } else if (*iordre == 2) {
1259         hermit[hermit_dim1 + 1] = .5;
1260         hermit[hermit_dim1 + 2] = -.9375;
1261         hermit[hermit_dim1 + 3] = 0.;
1262         hermit[hermit_dim1 + 4] = .625;
1263         hermit[hermit_dim1 + 5] = 0.;
1264         hermit[hermit_dim1 + 6] = -.1875;
1265
1266         hermit[(hermit_dim1 << 1) + 1] = .5;
1267         hermit[(hermit_dim1 << 1) + 2] = .9375;
1268         hermit[(hermit_dim1 << 1) + 3] = 0.;
1269         hermit[(hermit_dim1 << 1) + 4] = -.625;
1270         hermit[(hermit_dim1 << 1) + 5] = 0.;
1271         hermit[(hermit_dim1 << 1) + 6] = .1875;
1272
1273         hermit[hermit_dim1 * 3 + 1] = .3125;
1274         hermit[hermit_dim1 * 3 + 2] = -.4375;
1275         hermit[hermit_dim1 * 3 + 3] = -.375;
1276         hermit[hermit_dim1 * 3 + 4] = .625;
1277         hermit[hermit_dim1 * 3 + 5] = .0625;
1278         hermit[hermit_dim1 * 3 + 6] = -.1875;
1279
1280         hermit[(hermit_dim1 << 2) + 1] = -.3125;
1281         hermit[(hermit_dim1 << 2) + 2] = -.4375;
1282         hermit[(hermit_dim1 << 2) + 3] = .375;
1283         hermit[(hermit_dim1 << 2) + 4] = .625;
1284         hermit[(hermit_dim1 << 2) + 5] = -.0625;
1285         hermit[(hermit_dim1 << 2) + 6] = -.1875;
1286
1287         hermit[hermit_dim1 * 5 + 1] = .0625;
1288         hermit[hermit_dim1 * 5 + 2] = -.0625;
1289         hermit[hermit_dim1 * 5 + 3] = -.125;
1290         hermit[hermit_dim1 * 5 + 4] = .125;
1291         hermit[hermit_dim1 * 5 + 5] = .0625;
1292         hermit[hermit_dim1 * 5 + 6] = -.0625;
1293
1294         hermit[hermit_dim1 * 6 + 1] = .0625;
1295         hermit[hermit_dim1 * 6 + 2] = .0625;
1296         hermit[hermit_dim1 * 6 + 3] = -.125;
1297         hermit[hermit_dim1 * 6 + 4] = -.125;
1298         hermit[hermit_dim1 * 6 + 5] = .0625;
1299         hermit[hermit_dim1 * 6 + 6] = .0625;
1300     } else {
1301         *iercod = 1;
1302     }
1303
1304 /* ------------------------------ The End ------------------------------- 
1305 */
1306
1307     AdvApp2Var_SysBase::maermsg_("MMA1HER", iercod, 7L);
1308     if (ibb >= 3) {
1309         AdvApp2Var_SysBase::mgsomsg_("MMA1HER", 7L);
1310     }
1311     return 0;
1312 } /* mma1her_ */
1313 //=======================================================================
1314 //function : mma1jak_
1315 //purpose  : 
1316 //=======================================================================
1317 int mma1jak_(integer *ndimen, 
1318              integer *nbroot, 
1319              integer *iordre,
1320              integer *ndgjac, 
1321              doublereal *somtab, 
1322              doublereal *diftab, 
1323              doublereal *cgauss, 
1324              doublereal *crvjac, 
1325              integer *iercod)
1326 {
1327   /* System generated locals */
1328   integer somtab_dim1, somtab_offset, diftab_dim1, diftab_offset, 
1329   crvjac_dim1, crvjac_offset, cgauss_dim1;
1330
1331   /* Local variables */
1332   static integer ibb;
1333
1334 /* ********************************************************************** 
1335 */
1336
1337 /*     FUNCTION : */
1338 /*     ---------- */
1339 /*     Calculate the curve of approximation of a non-polynomial function */
1340 /*     in the base of Jacobi. */
1341
1342 /*     KEYWORDS : */
1343 /*     ----------- */
1344 /*     FUNCTION,DISCRETISATION,APPROXIMATION,CONSTRAINT,CURVE,JACOBI */
1345
1346 /*     INPUT ARGUMENTS : */
1347 /*     ------------------ */
1348 /*     NDIMEN: Total dimension of the space (sum of dimensions */
1349 /*             of sub-spaces) */
1350 /*     NBROOT: Nb of points of discretization of the iso, extremities not 
1351 /*             included. */
1352 /*     IORDRE: Order of constraint at the extremities of the boundary */
1353 /*              -1 = no constraints, */
1354 /*               0 = constraints of passage of limits (i.e. C0), */
1355 /*               1 = C0 + constraints of 1st derivatives (i.e. C1), */
1356 /*               2 = C1 + constraints of 2nd derivatives (i.e. C2). */
1357 /*     NDGJAC: Degree of development in series to be used for calculation in the  
1358 /*             base of Jacobi. */
1359
1360 /*     OUTPUT ARGUMENTS : */
1361 /*     ------------------- */
1362 /*     CRVJAC : Curve of approximation of FONCNP with (eventually) */
1363 /*              taking into account of constraints at the extremities. */
1364 /*              This curve is of degree NDGJAC. */
1365 /*     IERCOD : Error code : */
1366 /*               0 = All is ok. */
1367 /*              33 = Pb to return data of du block data */
1368 /*                   of coeff. of integration by GAUSS method. */
1369 /*                   by program MMAPPTT. */
1370
1371 /*     COMMONS USED   : */
1372 /*     ---------------- */
1373
1374 /*     REFERENCES CALLED   : */
1375 /*     ----------------------- */
1376 /* > */
1377 /* ********************************************************************** 
1378 */
1379
1380 /*   Name of the routine */
1381
1382     /* Parameter adjustments */
1383     diftab_dim1 = *nbroot / 2 + 1;
1384     diftab_offset = diftab_dim1;
1385     diftab -= diftab_offset;
1386     somtab_dim1 = *nbroot / 2 + 1;
1387     somtab_offset = somtab_dim1;
1388     somtab -= somtab_offset;
1389     crvjac_dim1 = *ndgjac + 1;
1390     crvjac_offset = crvjac_dim1;
1391     crvjac -= crvjac_offset;
1392     cgauss_dim1 = *nbroot / 2 + 1;
1393
1394     /* Function Body */
1395     ibb = AdvApp2Var_SysBase::mnfndeb_();
1396     if (ibb >= 2) {
1397         AdvApp2Var_SysBase::mgenmsg_("MMA1JAK", 7L);
1398     }
1399     *iercod = 0;
1400
1401 /* ----------------- Recover coeffs of integration by Gauss ----------- 
1402 */
1403
1404     AdvApp2Var_ApproxF2var::mmapptt_(ndgjac, nbroot, iordre, cgauss, iercod);
1405     if (*iercod > 0) {
1406         *iercod = 33;
1407         goto L9999;
1408     }
1409
1410 /* --------------- Calculate the curve in the base of Jacobi ----------- 
1411 */
1412
1413     mmmapcoe_(ndimen, ndgjac, iordre, nbroot, &somtab[somtab_offset], &diftab[
1414             diftab_offset], cgauss, &crvjac[crvjac_offset]);
1415
1416 /* ------------------------------ The End ------------------------------- 
1417 */
1418
1419 L9999:
1420     if (*iercod != 0) {
1421         AdvApp2Var_SysBase::maermsg_("MMA1JAK", iercod, 7L);
1422     }
1423     if (ibb >= 2) {
1424         AdvApp2Var_SysBase::mgsomsg_("MMA1JAK", 7L);
1425     }
1426     return 0;
1427 } /* mma1jak_ */
1428
1429 //=======================================================================
1430 //function : mma1noc_
1431 //purpose  : 
1432 //=======================================================================
1433 int mma1noc_(doublereal *dfuvin, 
1434              integer *ndimen, 
1435              integer *iordre, 
1436              doublereal *cntrin, 
1437              doublereal *duvout, 
1438              integer *isofav, 
1439              integer *ideriv, 
1440              doublereal *cntout)
1441 {
1442   /* System generated locals */
1443   integer i__1;
1444   doublereal d__1;
1445
1446   /* Local variables */
1447   static doublereal rider, riord;
1448   static integer nd, ibb;
1449   static doublereal bid;
1450 /* ********************************************************************** 
1451 */
1452
1453 /*     FUNCTION : */
1454 /*     ---------- */
1455 /*     Normalization of constraints of derivatives, defined on DFUVIN */
1456 /*     on block DUVOUT. */
1457
1458 /*     KEYWORDS : */
1459 /*     ----------- */
1460 /*     ALL, AB_SPECIFI::VECTEUR&,DERIVEE&,NORMALISATION,&VECTEUR */
1461
1462 /*     INPUT ARGUMENTS : */
1463 /*     ------------------ */
1464 /*     DFUVIN: Limits of the block of definition by U and by V where 
1465 */
1466 /*             constraints CNTRIN are defined. */
1467 /*     NDIMEN: Dimension of the space. */
1468 /*     IORDRE: Order of constraint imposed at the extremities of the iso. */
1469 /*             (if Iso-U, it is necessary to calculate derivatives by V and vice */
1470 /*             versa). */
1471 /*             = 0, the extremities of the iso are calculated */
1472 /*             = 1, additionally the 1st derivative in the direction */
1473 /*                  of the iso is calculated */
1474 /*             = 2, additionally the 2nd derivative in the direction */
1475 /*                  of the iso is calculated */
1476 /*     CNTRIN: Contains, if IORDRE>=0, IORDRE+1 derivatives */
1477 /*             of order IORDRE of F(Uc,v) or of F(u,Vc), following the */
1478 /*             value of ISOFAV, RENORMALIZED by u and v in (-1,1). */
1479 /*     DUVOUT: Limits of the block of definition by U and by V where the */
1480 /*             constraints CNTOUT will be defined. */
1481 /*     ISOFAV: Isoparameter fixed for the discretization; */
1482 /*             = 1, discretization with fixed U=Uc and variable V. */
1483 /*             = 2, discretization with fixed V=Vc and variable U. */
1484 /*     IDERIV: Ordre de derivee transverse a l'iso fixee (Si Iso-U=Uc */
1485 /*             is fixed, the derivative of order IDERIV is discretized by U */
1486 /*             of F(Uc,v). The same if iso-V is fixed). */
1487 /*             Varies from (positioning) to 2 (2nd derivative). */
1488
1489 /*     OUTPUT ARGUMENTS : */
1490 /*     ------------------- */
1491 /*     CNTOUT: Contains, if IORDRE>=0, IORDRE+1 derivatives */
1492 /*             of order IORDRE of F(Uc,v) or of F(u,Vc), depending on the */
1493 /*             value of ISOFAV, RENORMALIZED for u and v in DUVOUT. */
1494
1495 /*     COMMONS USED   : */
1496 /*     ---------------- */
1497
1498 /*     REFERENCES CALLED   : */
1499 /*     --------------------- */
1500
1501 /*     DESCRIPTION/NOTES/LIMITATIONS : */
1502 /*     ------------------------------- */
1503 /*     CNTRIN can be an output/input  argument, */
1504 /*     so the call: */
1505
1506 /*      CALL MMA1NOC(DFUVIN,NDIMEN,IORDRE,CNTRIN,DUVOUT */
1507 /*     1           ,ISOFAV,IDERIV,CNTRIN) */
1508
1509 /*     is correct. */
1510 /* > */
1511 /* ********************************************************************** 
1512 */
1513
1514 /*   Name of the routine */
1515
1516
1517     /* Parameter adjustments */
1518     dfuvin -= 3;
1519     --cntout;
1520     --cntrin;
1521     duvout -= 3;
1522
1523     /* Function Body */
1524     ibb = AdvApp2Var_SysBase::mnfndeb_();
1525     if (ibb >= 3) {
1526         AdvApp2Var_SysBase::mgenmsg_("MMA1NOC", 7L);
1527     }
1528
1529 /* --------------- Determination of coefficients of normalization -------
1530  */
1531
1532     if (*isofav == 1) {
1533         d__1 = (dfuvin[4] - dfuvin[3]) / (duvout[4] - duvout[3]);
1534         rider = AdvApp2Var_MathBase::pow__di(&d__1, ideriv);
1535         d__1 = (dfuvin[6] - dfuvin[5]) / (duvout[6] - duvout[5]);
1536         riord = AdvApp2Var_MathBase::pow__di(&d__1, iordre);
1537
1538     } else {
1539         d__1 = (dfuvin[6] - dfuvin[5]) / (duvout[6] - duvout[5]);
1540         rider = AdvApp2Var_MathBase::pow__di(&d__1, ideriv);
1541         d__1 = (dfuvin[4] - dfuvin[3]) / (duvout[4] - duvout[3]);
1542         riord = AdvApp2Var_MathBase::pow__di(&d__1, iordre);
1543     }
1544
1545 /* ------------- Renormalization of the vector of constraint --------------- 
1546 */
1547
1548     bid = rider * riord;
1549     i__1 = *ndimen;
1550     for (nd = 1; nd <= i__1; ++nd) {
1551         cntout[nd] = bid * cntrin[nd];
1552 /* L100: */
1553     }
1554
1555 /* ------------------------------ The end ------------------------------- 
1556 */
1557
1558     if (ibb >= 3) {
1559         AdvApp2Var_SysBase::mgsomsg_("MMA1NOC", 7L);
1560     }
1561     return 0;
1562 } /* mma1noc_ */
1563
1564 //=======================================================================
1565 //function : mma1nop_
1566 //purpose  : 
1567 //=======================================================================
1568 int mma1nop_(integer *nbroot, 
1569              doublereal *rootlg, 
1570              doublereal *uvfonc, 
1571              integer *isofav, 
1572              doublereal *ttable, 
1573              integer *iercod)
1574
1575 {
1576   /* System generated locals */
1577   integer i__1;
1578
1579   /* Local variables */
1580   static doublereal alinu, blinu, alinv, blinv;
1581   static integer ii, ibb;
1582
1583 /* ***********************************************************************
1584  */
1585
1586 /*     FUNCTION : */
1587 /*     ---------- */
1588 /*     Normalization of parameters of an iso, starting from  */
1589 /*     parametric block and parameters on (-1,1). */
1590
1591 /*     KEYWORDS : */
1592 /*     ----------- */
1593 /*      TOUS,AB_SPECIFI :: ISO&,POINT&,NORMALISATION,&POINT,&ISO */
1594
1595 /*     INPUT ARGUMENTS : */
1596 /*     -------------------- */
1597 /*        NBROOT: Nb of points of discretisation INSIDE the iso */
1598 /*                defined on (-1,1). */
1599 /*        ROOTLG: Table of discretization parameters on )-1,1( */
1600 /*                of the iso. */
1601 /*        UVFONC: Block of definition of the iso */
1602 /*        ISOFAV: = 1, this is iso-u; =2, this is iso-v. */
1603
1604 /*     OUTPUT ARGUMENTS : */
1605 /*     --------------------- */
1606 /*        TTABLE: Table of parameters renormalized on UVFONC of the iso. 
1607 */
1608 /*        IERCOD: = 0, OK */
1609 /*                = 1, ISOFAV is out of allowed values. */
1610
1611 /* > */
1612 /* ********************************************************************** 
1613 */
1614 /*   Name of the routine */
1615
1616
1617     /* Parameter adjustments */
1618     --rootlg;
1619     uvfonc -= 3;
1620
1621     /* Function Body */
1622     ibb = AdvApp2Var_SysBase::mnfndeb_();
1623     if (ibb >= 3) {
1624         AdvApp2Var_SysBase::mgenmsg_("MMA1NOP", 7L);
1625     }
1626
1627     alinu = (uvfonc[4] - uvfonc[3]) / 2.;
1628     blinu = (uvfonc[4] + uvfonc[3]) / 2.;
1629     alinv = (uvfonc[6] - uvfonc[5]) / 2.;
1630     blinv = (uvfonc[6] + uvfonc[5]) / 2.;
1631
1632     if (*isofav == 1) {
1633         ttable[0] = uvfonc[5];
1634         i__1 = *nbroot;
1635         for (ii = 1; ii <= i__1; ++ii) {
1636             ttable[ii] = alinv * rootlg[ii] + blinv;
1637 /* L100: */
1638         }
1639         ttable[*nbroot + 1] = uvfonc[6];
1640     } else if (*isofav == 2) {
1641         ttable[0] = uvfonc[3];
1642         i__1 = *nbroot;
1643         for (ii = 1; ii <= i__1; ++ii) {
1644             ttable[ii] = alinu * rootlg[ii] + blinu;
1645 /* L200: */
1646         }
1647         ttable[*nbroot + 1] = uvfonc[4];
1648     } else {
1649         goto L9100;
1650     }
1651
1652     goto L9999;
1653
1654 /* ------------------------------ THE END ------------------------------- 
1655 */
1656
1657 L9100:
1658     *iercod = 1;
1659     goto L9999;
1660
1661 L9999:
1662     if (*iercod != 0) {
1663         AdvApp2Var_SysBase::maermsg_("MMA1NOP", iercod, 7L);
1664     }
1665     if (ibb >= 3) {
1666         AdvApp2Var_SysBase::mgsomsg_("MMA1NOP", 7L);
1667     }
1668
1669  return 0 ;
1670
1671 } /* mma1nop_ */
1672
1673 //=======================================================================
1674 //function : mma2ac1_
1675 //purpose  : 
1676 //=======================================================================
1677 int AdvApp2Var_ApproxF2var::mma2ac1_(integer const *ndimen, 
1678                                      integer const *mxujac, 
1679                                      integer const *mxvjac, 
1680                                      integer const *iordru, 
1681                                      integer const *iordrv, 
1682                                      doublereal const *contr1, 
1683                                      doublereal const * contr2, 
1684                                      doublereal const *contr3, 
1685                                      doublereal const *contr4, 
1686                                      doublereal const *uhermt, 
1687                                      doublereal const *vhermt, 
1688                                      doublereal *patjac)
1689
1690 {
1691   /* System generated locals */
1692   integer contr1_dim1, contr1_dim2, contr1_offset, contr2_dim1, contr2_dim2,
1693   contr2_offset, contr3_dim1, contr3_dim2, contr3_offset, 
1694   contr4_dim1, contr4_dim2, contr4_offset, uhermt_dim1, 
1695   uhermt_offset, vhermt_dim1, vhermt_offset, patjac_dim1, 
1696   patjac_dim2, patjac_offset, i__1, i__2, i__3, i__4, i__5;
1697
1698   /* Local variables */
1699   static logical ldbg;
1700   static integer ndgu, ndgv;
1701   static doublereal bidu1, bidu2, bidv1, bidv2;
1702   static integer ioru1, iorv1, ii, nd, jj, ku, kv;
1703   static doublereal cnt1, cnt2, cnt3, cnt4;
1704
1705 /* ********************************************************************** 
1706 */
1707
1708 /*     FUNCTION : */
1709 /*     ---------- */
1710 /*     Add polynoms of edge constraints. */
1711
1712 /*     KEYWORDS : */
1713 /*     ----------- */
1714 /*  TOUS,AB_SPECIFI::POINT&,CONTRAINTE&,ADDITION,&POLYNOME */
1715
1716 /*     INPUT ARGUMENTS  : */
1717 /*     ------------------ */
1718 /*   NDIMEN: Dimension of the space. */
1719 /*   MXUJAC: Max degree of the polynom of approximation by U. The  */
1720 /*           representation in the orthogonal base starts from degree */
1721 /*           0 to degree MXUJAC-2*(IORDRU+1). The polynomial base is the */
1722 /*           base of Jacobi of order -1 (Legendre), 0, 1 or 2 */
1723 /*   MXVJAC: Max degree of the polynom of approximation by V. 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 /*   IORDRU: Order of the base of Jacobi (-1,0,1 or 2) by U. Corresponds */
1728 /*           to the step of constraints: C0, C1 or C2. */
1729 /*   IORDRV: Order of the base of Jacobi (-1,0,1 or 2) by V. Corresponds */
1730 /*           to the step of constraints: C0, C1 or C2. */
1731 /*   CONTR1: Contains, if IORDRU and IORDRV>=0, the values at the */
1732 /*           extremities of F(U0,V0) and its derivatives. */
1733 /*   CONTR2: Contains, if IORDRU and IORDRV>=0, the values at the */
1734 /*           extremities of F(U1,V0) and its derivatives. */
1735 /*   CONTR3: Contains, if IORDRU and IORDRV>=0, the values at the */
1736 /*           extremities of F(U0,V1) and its derivatives. */
1737 /*   CONTR4: Contains, if IORDRU and IORDRV>=0, the values at the */
1738 /*           extremities of F(U1,V1) and its derivatives. */
1739 /*   UHERMT: Coeff. of Hermit polynoms of order IORDRU. */
1740 /*   VHERMT: Coeff. of Hermit polynoms of order IORDRV. */
1741 /*   PATJAC: Table of coefficients of the polynom P(u,v) of approximation */
1742 /*           of F(u,v) WITHOUT taking into account the constraints. */
1743
1744 /*     OUTPUT ARGUMENTS : */
1745 /*     ------------------- */
1746 /*   PATJAC: Table of coefficients of the polynom P(u,v) by approximation */
1747 /*           of F(u,v) WITH taking into account of constraints. */
1748 /* > */
1749 /* ********************************************************************** 
1750 */
1751 /*   Name of the routine */
1752
1753 /* --------------------------- Initialization -------------------------- 
1754 */
1755
1756     /* Parameter adjustments */
1757     patjac_dim1 = *mxujac + 1;
1758     patjac_dim2 = *mxvjac + 1;
1759     patjac_offset = patjac_dim1 * patjac_dim2;
1760     patjac -= patjac_offset;
1761     uhermt_dim1 = (*iordru << 1) + 2;
1762     uhermt_offset = uhermt_dim1;
1763     uhermt -= uhermt_offset;
1764     vhermt_dim1 = (*iordrv << 1) + 2;
1765     vhermt_offset = vhermt_dim1;
1766     vhermt -= vhermt_offset;
1767     contr4_dim1 = *ndimen;
1768     contr4_dim2 = *iordru + 2;
1769     contr4_offset = contr4_dim1 * (contr4_dim2 + 1) + 1;
1770     contr4 -= contr4_offset;
1771     contr3_dim1 = *ndimen;
1772     contr3_dim2 = *iordru + 2;
1773     contr3_offset = contr3_dim1 * (contr3_dim2 + 1) + 1;
1774     contr3 -= contr3_offset;
1775     contr2_dim1 = *ndimen;
1776     contr2_dim2 = *iordru + 2;
1777     contr2_offset = contr2_dim1 * (contr2_dim2 + 1) + 1;
1778     contr2 -= contr2_offset;
1779     contr1_dim1 = *ndimen;
1780     contr1_dim2 = *iordru + 2;
1781     contr1_offset = contr1_dim1 * (contr1_dim2 + 1) + 1;
1782     contr1 -= contr1_offset;
1783
1784     /* Function Body */
1785     ldbg = AdvApp2Var_SysBase::mnfndeb_() >= 3;
1786     if (ldbg) {
1787         AdvApp2Var_SysBase::mgenmsg_("MMA2AC1", 7L);
1788     }
1789
1790 /* ------------ SUBTRACTION OF ANGULAR CONSTRAINTS ------------------- 
1791 */
1792
1793     ioru1 = *iordru + 1;
1794     iorv1 = *iordrv + 1;
1795     ndgu = (*iordru << 1) + 1;
1796     ndgv = (*iordrv << 1) + 1;
1797
1798     i__1 = iorv1;
1799     for (jj = 1; jj <= i__1; ++jj) {
1800         i__2 = ioru1;
1801         for (ii = 1; ii <= i__2; ++ii) {
1802             i__3 = *ndimen;
1803             for (nd = 1; nd <= i__3; ++nd) {
1804                 cnt1 = contr1[nd + (ii + jj * contr1_dim2) * contr1_dim1];
1805                 cnt2 = contr2[nd + (ii + jj * contr2_dim2) * contr2_dim1];
1806                 cnt3 = contr3[nd + (ii + jj * contr3_dim2) * contr3_dim1];
1807                 cnt4 = contr4[nd + (ii + jj * contr4_dim2) * contr4_dim1];
1808                 i__4 = ndgv;
1809                 for (kv = 0; kv <= i__4; ++kv) {
1810                     bidv1 = vhermt[kv + ((jj << 1) - 1) * vhermt_dim1];
1811                     bidv2 = vhermt[kv + (jj << 1) * vhermt_dim1];
1812                     i__5 = ndgu;
1813                     for (ku = 0; ku <= i__5; ++ku) {
1814                         bidu1 = uhermt[ku + ((ii << 1) - 1) * uhermt_dim1];
1815                         bidu2 = uhermt[ku + (ii << 1) * uhermt_dim1];
1816                         patjac[ku + (kv + nd * patjac_dim2) * patjac_dim1] = 
1817                                 patjac[ku + (kv + nd * patjac_dim2) * 
1818                                 patjac_dim1] - bidu1 * bidv1 * cnt1 - bidu2 * 
1819                                 bidv1 * cnt2 - bidu1 * bidv2 * cnt3 - bidu2 * 
1820                                 bidv2 * cnt4;
1821 /* L500: */
1822                     }
1823 /* L400: */
1824                 }
1825 /* L300: */
1826             }
1827 /* L200: */
1828         }
1829 /* L100: */
1830     }
1831
1832 /* ------------------------------ The end ------------------------------- 
1833 */
1834
1835     if (ldbg) {
1836         AdvApp2Var_SysBase::mgsomsg_("MMA2AC1", 7L);
1837     }
1838     return 0;
1839 } /* mma2ac1_ */
1840
1841 //=======================================================================
1842 //function : mma2ac2_
1843 //purpose  : 
1844 //=======================================================================
1845 int AdvApp2Var_ApproxF2var::mma2ac2_(const integer *ndimen, 
1846                                      const integer *mxujac, 
1847                                      const integer *mxvjac, 
1848                                      const integer *iordrv, 
1849                                      const integer *nclimu, 
1850                                      const integer *ncfiv1, 
1851                                      const doublereal *crbiv1, 
1852                                      const integer *ncfiv2, 
1853                                      const doublereal *crbiv2, 
1854                                      const doublereal *vhermt, 
1855                                      doublereal *patjac)
1856
1857 {
1858   /* System generated locals */
1859   integer crbiv1_dim1, crbiv1_dim2, crbiv1_offset, crbiv2_dim1, crbiv2_dim2,
1860   crbiv2_offset, patjac_dim1, patjac_dim2, patjac_offset, 
1861   vhermt_dim1, vhermt_offset, i__1, i__2, i__3, i__4;
1862
1863   /* Local variables */
1864   static logical ldbg;
1865   static integer ndgv1, ndgv2, ii, jj, nd, kk;
1866   static doublereal bid1, bid2;
1867
1868 /* ********************************************************************** 
1869 */
1870
1871 /*     FUNCTION : */
1872 /*     ---------- */
1873 /*     Add polynoms of constraints */
1874
1875 /*     KEYWORDS : */
1876 /*     ----------- */
1877 /*     FUNCTION,APPROXIMATION,COEFFICIENT,POLYNOM */
1878
1879 /*     INPUT ARGUMENTS : */
1880 /*     ------------------ */
1881 /*   NDIMEN: Dimension of the space. */
1882 /*   MXUJAC: Max degree of the polynom of approximation by U. The  */
1883 /*           representation in the orthogonal base starts from degree */
1884 /*           0 to degree MXUJAC-2*(IORDRU+1). The polynomial base is the */
1885 /*           base of Jacobi of order -1 (Legendre), 0, 1 or 2 */
1886 /*   MXVJAC: Max degree of the polynom of approximation by V. 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 /*   IORDRV: Order of the base of Jacobi (-1,0,1 or 2) by V. Corresponds */
1891 /*           to the step of constraints: C0, C1 or C2. */
1892 /*   NCLIMU  LIMIT nb of coeff by u of the solution P(u,v) 
1893 *    NCFIV1: Nb of Coeff. of curves stored in CRBIV1. */
1894 /*   CRBIV1: Table of coeffs of the approximation of iso-V0 and its */
1895 /*           derivatives till order IORDRV. */
1896 /*   NCFIV2: Nb of Coeff. of curves stored in CRBIV2. */
1897 /*   CRBIV2: Table of coeffs of approximation of iso-V1 and its */
1898 /*           derivatives till order IORDRV. */
1899 /*   VHERMT: Coeff. of Hermit polynoms of order IORDRV. */
1900 /*   PATJAC: Table of coefficients of the polynom P(u,v) of approximation */
1901 /*           of F(u,v) WITHOUT taking into account the constraints. */
1902
1903 /*     OUTPUT ARGUMENTS : */
1904 /*     ------------------- */
1905 /*   PATJAC: Table of coefficients of the polynom P(u,v) by approximation */
1906 /*           of F(u,v) WITH taking into account of constraints. */
1907 /* > *//*
1908
1909
1910 /* > */
1911 /* ********************************************************************** 
1912 */
1913 /*   Name of the routine */
1914
1915 /* --------------------------- Initialisations -------------------------- 
1916 */
1917
1918     /* Parameter adjustments */
1919     patjac_dim1 = *mxujac + 1;
1920     patjac_dim2 = *mxvjac + 1;
1921     patjac_offset = patjac_dim1 * patjac_dim2;
1922     patjac -= patjac_offset;
1923     vhermt_dim1 = (*iordrv << 1) + 2;
1924     vhermt_offset = vhermt_dim1;
1925     vhermt -= vhermt_offset;
1926     --ncfiv2;
1927     --ncfiv1;
1928     crbiv2_dim1 = *nclimu;
1929     crbiv2_dim2 = *ndimen;
1930     crbiv2_offset = crbiv2_dim1 * (crbiv2_dim2 + 1);
1931     crbiv2 -= crbiv2_offset;
1932     crbiv1_dim1 = *nclimu;
1933     crbiv1_dim2 = *ndimen;
1934     crbiv1_offset = crbiv1_dim1 * (crbiv1_dim2 + 1);
1935     crbiv1 -= crbiv1_offset;
1936
1937     /* Function Body */
1938     ldbg = AdvApp2Var_SysBase::mnfndeb_() >= 3;
1939     if (ldbg) {
1940         AdvApp2Var_SysBase::mgenmsg_("MMA2AC2", 7L);
1941     }
1942
1943 /* ------------ ADDING of coeff by u of curves, by v of Hermit -------- 
1944 */
1945
1946     i__1 = *iordrv + 1;
1947     for (ii = 1; ii <= i__1; ++ii) {
1948         ndgv1 = ncfiv1[ii] - 1;
1949         ndgv2 = ncfiv2[ii] - 1;
1950         i__2 = *ndimen;
1951         for (nd = 1; nd <= i__2; ++nd) {
1952             i__3 = (*iordrv << 1) + 1;
1953             for (jj = 0; jj <= i__3; ++jj) {
1954                 bid1 = vhermt[jj + ((ii << 1) - 1) * vhermt_dim1];
1955                 i__4 = ndgv1;
1956                 for (kk = 0; kk <= i__4; ++kk) {
1957                     patjac[kk + (jj + nd * patjac_dim2) * patjac_dim1] += 
1958                             bid1 * crbiv1[kk + (nd + ii * crbiv1_dim2) * 
1959                             crbiv1_dim1];
1960 /* L400: */
1961                 }
1962                 bid2 = vhermt[jj + (ii << 1) * vhermt_dim1];
1963                 i__4 = ndgv2;
1964                 for (kk = 0; kk <= i__4; ++kk) {
1965                     patjac[kk + (jj + nd * patjac_dim2) * patjac_dim1] += 
1966                             bid2 * crbiv2[kk + (nd + ii * crbiv2_dim2) * 
1967                             crbiv2_dim1];
1968 /* L500: */
1969                 }
1970 /* L300: */
1971             }
1972 /* L200: */
1973         }
1974 /* L100: */
1975     }
1976
1977 /* ------------------------------ The end ------------------------------- 
1978 */
1979
1980     if (ldbg) {
1981         AdvApp2Var_SysBase::mgsomsg_("MMA2AC2", 7L);
1982     }
1983     return 0;
1984 } /* mma2ac2_ */
1985
1986
1987 //=======================================================================
1988 //function : mma2ac3_
1989 //purpose  : 
1990 //=======================================================================
1991 int AdvApp2Var_ApproxF2var::mma2ac3_(const integer *ndimen, 
1992                                      const integer *mxujac, 
1993                                      const integer *mxvjac, 
1994                                      const integer *iordru, 
1995                                      const integer *nclimv, 
1996                                      const integer *ncfiu1, 
1997                                      const doublereal * crbiu1, 
1998                                      const integer *ncfiu2, 
1999                                      const doublereal *crbiu2, 
2000                                      const doublereal *uhermt, 
2001                                      doublereal *patjac)
2002
2003 {
2004   /* System generated locals */
2005   integer crbiu1_dim1, crbiu1_dim2, crbiu1_offset, crbiu2_dim1, crbiu2_dim2,
2006   crbiu2_offset, patjac_dim1, patjac_dim2, patjac_offset, 
2007   uhermt_dim1, uhermt_offset, i__1, i__2, i__3, i__4;
2008
2009   /* Local variables */
2010   static logical ldbg;
2011   static integer ndgu1, ndgu2, ii, jj, nd, kk;
2012   static doublereal bid1, bid2;
2013
2014 /* ********************************************************************** 
2015 */
2016
2017 /*     FUNCTION : */
2018 /*     ---------- */
2019 /*     Ajout des polynomes de contraintes */
2020
2021 /*     KEYWORDS : */
2022 /*     ----------- */
2023 /*     FONCTION,APPROXIMATION,COEFFICIENT,POLYNOME */
2024
2025 /*     INPUT ARGUMENTS : */
2026 /*     ------------------ */
2027 /*   NDIMEN: Dimension of the space. */
2028 /*   MXUJAC: Max degree of the polynom of approximation by U. The  */
2029 /*           representation in the orthogonal base starts from degree */
2030 /*           0 to degree MXUJAC-2*(IORDRU+1). The polynomial base is the */
2031 /*           base of Jacobi of order -1 (Legendre), 0, 1 or 2 */
2032 /*   MXVJAC: Max degree of the polynom of approximation by V. 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 /*   IORDRU: Order of the base of Jacobi (-1,0,1 or 2) by U. Corresponds */
2037 /*           to the step of constraints: C0, C1 or C2. */
2038 /*   NCLIMV  LIMIT nb of coeff by v of the solution P(u,v) 
2039 *    NCFIU1: Nb of Coeff. of curves stored in CRBIU1. */
2040 /*   CRBIU1: Table of coeffs of the approximation of iso-U0 and its */
2041 /*           derivatives till order IORDRU. */
2042 /*   NCFIU2: Nb of Coeff. of curves stored in CRBIU2. */
2043 /*   CRBIU2: Table of coeffs of approximation of iso-U1 and its */
2044 /*           derivatives till order IORDRU */
2045 /*   UHERMT: Coeff. of Hermit polynoms of order IORDRU. */
2046 /*   PATJAC: Table of coefficients of the polynom P(u,v) of approximation */
2047 /*           of F(u,v) WITHOUT taking into account the constraints. */
2048
2049 /*     OUTPUT ARGUMENTS : */
2050 /*     ------------------- */
2051 /*   PATJAC: Table of coefficients of the polynom P(u,v) by approximation */
2052 /*           of F(u,v) WITH taking into account of constraints. */
2053
2054
2055 /* > */
2056 /* ********************************************************************** 
2057 */
2058 /*   The name of the routine */
2059
2060 /* --------------------------- Initializations -------------------------- 
2061 */
2062
2063     /* Parameter adjustments */
2064     patjac_dim1 = *mxujac + 1;
2065     patjac_dim2 = *mxvjac + 1;
2066     patjac_offset = patjac_dim1 * patjac_dim2;
2067     patjac -= patjac_offset;
2068     uhermt_dim1 = (*iordru << 1) + 2;
2069     uhermt_offset = uhermt_dim1;
2070     uhermt -= uhermt_offset;
2071     --ncfiu2;
2072     --ncfiu1;
2073     crbiu2_dim1 = *nclimv;
2074     crbiu2_dim2 = *ndimen;
2075     crbiu2_offset = crbiu2_dim1 * (crbiu2_dim2 + 1);
2076     crbiu2 -= crbiu2_offset;
2077     crbiu1_dim1 = *nclimv;
2078     crbiu1_dim2 = *ndimen;
2079     crbiu1_offset = crbiu1_dim1 * (crbiu1_dim2 + 1);
2080     crbiu1 -= crbiu1_offset;
2081
2082     /* Function Body */
2083     ldbg = AdvApp2Var_SysBase::mnfndeb_() >= 3;
2084     if (ldbg) {
2085         AdvApp2Var_SysBase::mgenmsg_("MMA2AC3", 7L);
2086     }
2087
2088 /* ------------ ADDING of coeff by u of curves, by v of Hermit -------- 
2089 */
2090
2091     i__1 = *iordru + 1;
2092     for (ii = 1; ii <= i__1; ++ii) {
2093         ndgu1 = ncfiu1[ii] - 1;
2094         ndgu2 = ncfiu2[ii] - 1;
2095         i__2 = *ndimen;
2096         for (nd = 1; nd <= i__2; ++nd) {
2097             i__3 = ndgu1;
2098             for (jj = 0; jj <= i__3; ++jj) {
2099                 bid1 = crbiu1[jj + (nd + ii * crbiu1_dim2) * crbiu1_dim1];
2100                 i__4 = (*iordru << 1) + 1;
2101                 for (kk = 0; kk <= i__4; ++kk) {
2102                     patjac[kk + (jj + nd * patjac_dim2) * patjac_dim1] += 
2103                             bid1 * uhermt[kk + ((ii << 1) - 1) * uhermt_dim1];
2104 /* L400: */
2105                 }
2106 /* L300: */
2107             }
2108             i__3 = ndgu2;
2109             for (jj = 0; jj <= i__3; ++jj) {
2110                 bid2 = crbiu2[jj + (nd + ii * crbiu2_dim2) * crbiu2_dim1];
2111                 i__4 = (*iordru << 1) + 1;
2112                 for (kk = 0; kk <= i__4; ++kk) {
2113                     patjac[kk + (jj + nd * patjac_dim2) * patjac_dim1] += 
2114                             bid2 * uhermt[kk + (ii << 1) * uhermt_dim1];
2115 /* L600: */
2116                 }
2117 /* L500: */
2118             }
2119
2120 /* L200: */
2121         }
2122 /* L100: */
2123     }
2124
2125 /* ------------------------------ The end ------------------------------- 
2126 */
2127
2128     if (ldbg) {
2129         AdvApp2Var_SysBase::mgsomsg_("MMA2AC3", 7L);
2130     }
2131     return 0;
2132 } /* mma2ac3_ */
2133
2134 //=======================================================================
2135 //function : mma2can_
2136 //purpose  : 
2137 //=======================================================================
2138 int AdvApp2Var_ApproxF2var::mma2can_(const integer *ncfmxu, 
2139                                      const integer *ncfmxv,
2140                                      const integer *ndimen, 
2141                                      const integer *iordru, 
2142                                      const integer *iordrv, 
2143                                      const integer *ncoefu, 
2144                                      const integer *ncoefv, 
2145                                      const doublereal *patjac, 
2146                                      doublereal *pataux, 
2147                                      doublereal *patcan, 
2148                                      integer *iercod)
2149
2150 {
2151   /* System generated locals */
2152   integer patjac_dim1, patjac_dim2, patjac_offset, patcan_dim1, patcan_dim2,
2153   patcan_offset, i__1, i__2;
2154
2155   /* Local variables */
2156   static logical ldbg;
2157   static integer ilon1, ilon2, ii, nd;
2158
2159 /* ********************************************************************** 
2160 */
2161
2162 /*     FUNCTION : */
2163 /*     ---------- */
2164 /*     Change of Jacobi base to canonical (-1,1) and writing in a greater */
2165 /*     table. */
2166
2167 /*     KEYWORDS : */
2168 /*     ----------- */
2169 /*     ALL,AB_SPECIFI,CARREAU&,CONVERSION,JACOBI,CANNONIQUE,&CARREAU */
2170
2171 /*     INPUT ARGUMENTS : */
2172 /*     -------------------- */
2173 /*     NCFMXU: Dimension by U of resulting table PATCAN */
2174 /*     NCFMXV: Dimension by V of resulting table PATCAN */
2175 /*     NDIMEN: Dimension of the workspace. */
2176 /*     IORDRU: Order of constraint by U */
2177 /*     IORDRV: Order of constraint by V. */
2178 /*     NCOEFU: Nb of coeff by U of square PATJAC */
2179 /*     NCOEFV: Nb of coeff by V of square PATJAC */
2180 /*     PATJAC: Square in the base of Jacobi of order IORDRU by U and */
2181 /*             IORDRV by V. */
2182
2183 /*     OUTPUT ARGUMENTS : */
2184 /*     --------------------- */
2185 /*     PATAUX: Auxiliary Table. */
2186 /*     PATCAN: Table of coefficients in the canonic base. */
2187 /*     IERCOD: Error code. */
2188 /*             = 0, everything goes well, and all things are equal. */
2189 /*             = 1, the program refuses to process with incorrect input arguments */
2190
2191
2192 /*     COMMONS USED : */
2193 /*     ------------------ */
2194
2195 /*     REFERENCES CALLED : */
2196 /*     --------------------- */
2197
2198 /*     DESCRIPTION/NOTES/LIMITATIONS : */
2199 /*     ----------------------------------- */
2200 /* > */
2201 /* ********************************************************************** 
2202 */
2203
2204
2205     /* Parameter adjustments */
2206     patcan_dim1 = *ncfmxu;
2207     patcan_dim2 = *ncfmxv;
2208     patcan_offset = patcan_dim1 * (patcan_dim2 + 1) + 1;
2209     patcan -= patcan_offset;
2210     --pataux;
2211     patjac_dim1 = *ncoefu;
2212     patjac_dim2 = *ncoefv;
2213     patjac_offset = patjac_dim1 * (patjac_dim2 + 1) + 1;
2214     patjac -= patjac_offset;
2215
2216     /* Function Body */
2217     ldbg = AdvApp2Var_SysBase::mnfndeb_() >= 2;
2218     if (ldbg) {
2219         AdvApp2Var_SysBase::mgenmsg_("MMA2CAN", 7L);
2220     }
2221     *iercod = 0;
2222
2223     if (*iordru < -1 || *iordru > 2) {
2224         goto L9100;
2225     }
2226     if (*iordrv < -1 || *iordrv > 2) {
2227         goto L9100;
2228     }
2229     if (*ncoefu > *ncfmxu || *ncoefv > *ncfmxv) {
2230         goto L9100;
2231     }
2232
2233 /* --> Pass to canonic base (-1,1) */
2234     mmjacpt_(ndimen, ncoefu, ncoefv, iordru, iordrv, &patjac[patjac_offset], &
2235             pataux[1], &patcan[patcan_offset]);
2236
2237 /* --> Write all in a greater table */
2238     AdvApp2Var_MathBase::mmfmca8_((integer *)ncoefu, 
2239              (integer *)ncoefv, 
2240              (integer *)ndimen, 
2241              (integer *)ncfmxu, 
2242              (integer *)ncfmxv, 
2243              (integer *)ndimen, 
2244              (doublereal *)&patcan[patcan_offset], 
2245              (doublereal *)&patcan[patcan_offset]);
2246
2247 /* --> Complete with zeros the resulting table. */
2248     ilon1 = *ncfmxu - *ncoefu;
2249     ilon2 = *ncfmxu * (*ncfmxv - *ncoefv);
2250     i__1 = *ndimen;
2251     for (nd = 1; nd <= i__1; ++nd) {
2252         if (ilon1 > 0) {
2253             i__2 = *ncoefv;
2254             for (ii = 1; ii <= i__2; ++ii) {
2255                 AdvApp2Var_SysBase::mvriraz_(&ilon1, 
2256                          (char *)&patcan[*ncoefu + 1 + (ii + nd * patcan_dim2) * patcan_dim1]);
2257 /* L110: */
2258             }
2259         }
2260         if (ilon2 > 0) {
2261             AdvApp2Var_SysBase::mvriraz_(&ilon2, 
2262                      (char *)&patcan[(*ncoefv + 1 + nd * patcan_dim2) * patcan_dim1 + 1]);
2263         }
2264 /* L100: */
2265     }
2266
2267     goto L9999;
2268
2269 /* ---------------------- 
2270 */
2271
2272 L9100:
2273     *iercod = 1;
2274     goto L9999;
2275
2276 L9999:
2277     AdvApp2Var_SysBase::maermsg_("MMA2CAN", iercod, 7L);
2278     if (ldbg) {
2279         AdvApp2Var_SysBase::mgsomsg_("MMA2CAN", 7L);
2280     }
2281  return 0 ;
2282 } /* mma2can_ */
2283
2284 //=======================================================================
2285 //function : mma2cd1_
2286 //purpose  : 
2287 //=======================================================================
2288 int mma2cd1_(integer *ndimen, 
2289              integer *nbpntu, 
2290              doublereal *urootl, 
2291              integer *nbpntv, 
2292              doublereal *vrootl, 
2293              integer *iordru, 
2294              integer *iordrv, 
2295              doublereal *contr1, 
2296              doublereal *contr2, 
2297              doublereal *contr3, 
2298              doublereal *contr4, 
2299              doublereal *fpntbu, 
2300              doublereal *fpntbv, 
2301              doublereal *uhermt, 
2302              doublereal *vhermt, 
2303              doublereal *sosotb, 
2304              doublereal *soditb, 
2305              doublereal *disotb, 
2306              doublereal *diditb)
2307
2308 {
2309   static integer c__1 = 1;
2310
2311 /* System generated locals */
2312     integer contr1_dim1, contr1_dim2, contr1_offset, contr2_dim1, contr2_dim2,
2313              contr2_offset, contr3_dim1, contr3_dim2, contr3_offset, 
2314             contr4_dim1, contr4_dim2, contr4_offset, uhermt_dim1, 
2315             uhermt_offset, vhermt_dim1, vhermt_offset, fpntbu_dim1, 
2316             fpntbu_offset, fpntbv_dim1, fpntbv_offset, sosotb_dim1, 
2317             sosotb_dim2, sosotb_offset, diditb_dim1, diditb_dim2, 
2318             diditb_offset, soditb_dim1, soditb_dim2, soditb_offset, 
2319             disotb_dim1, disotb_dim2, disotb_offset, i__1, i__2, i__3, i__4, 
2320             i__5;
2321
2322     /* Local variables */
2323     static integer ncfhu, ncfhv, nuroo, nvroo, nd, ii, jj, kk, ll, ibb, kkm, 
2324             llm, kkp, llp;
2325     static doublereal bid1, bid2, bid3, bid4;
2326     static doublereal diu1, diu2, div1, div2, sou1, sou2, sov1, sov2;
2327
2328 /* ********************************************************************** 
2329 */
2330
2331 /*     FUNCTION : */
2332 /*     ---------- */
2333 /*     Discretisation on the parameters of polynoms of interpolation */
2334 /*     of constraints at the corners of order IORDRE. */
2335
2336 /*     KEYWORDS : */
2337 /*     ----------- */
2338 /*     TOUS, AB_SPECIFI::CONTRAINTE&, DISCRETISATION, &POINT */
2339
2340 /*     INPUT ARGUMENTS : */
2341 /*     ------------------ */
2342 /*     NDIMEN: Dimension of the space. */
2343 /*     NBPNTU: Nb of INTERNAL parameters of discretisation by U. */
2344 /*             This is also the nb of root of Legendre polynom where discretization is done. */
2345 /*     UROOTL: Table of parameters of discretisation ON (-1,1) by U. 
2346 */
2347 /*     NBPNTV: Nb of INTERNAL  parameters of discretisation by V. */
2348 /*             This is also the nb of root of Legendre polynom where discretization is done. */
2349 /*     VROOTL: Table of discretization parameters on (-1,1) by V. 
2350 /*     IORDRU: Order of constraint imposed at the extremities of iso-V */
2351 /*             = 0, calculate the extremities of iso-V */
2352 /*             = 1, calculate, additionally, the 1st derivative in the direction of iso-V */
2353 /*             = 2, calculate, additionally, the 2nd derivative in the direction of iso-V */
2354 /*     IORDRV: Order of constraint imposed at the extremities of iso-U */
2355 /*             = 0, calculate the extremities of iso-U */
2356 /*             = 1, calculate, additionally, the 1st derivative in the direction of iso-U */
2357 /*             = 2, calculate, additionally, the 2nd derivative in the direction of iso-U */
2358 /*   CONTR1: Contains, if IORDRU and IORDRV>=0, the values at the */
2359 /*           extremities of F(U0,V0) and its derivatives. */
2360 /*   CONTR2: Contains, if IORDRU and IORDRV>=0, the values at the */
2361 /*           extremities of F(U1,V0) and its derivatives. */
2362 /*   CONTR3: Contains, if IORDRU and IORDRV>=0, the values at the */
2363 /*           extremities of F(U0,V1) and its derivatives. */
2364 /*   CONTR4: Contains, if IORDRU and IORDRV>=0, the values at the */
2365 /*           extremities of F(U1,V1) and its derivatives. */
2366 /*     SOSOTB: Preinitialized table (input/output argument). */
2367 /*     DISOTB: Preinitialized table (input/output argument). */
2368 /*     SODITB: Preinitialized table (input/output argument). */
2369 /*     DIDITB: Preinitialized table (input/output argument) */
2370
2371 /*     OUTPUT ARGUMENTS : */
2372 /*     ------------------- */
2373 /*     FPNTBU: Auxiliary table. */
2374 /*     FPNTBV: Auxiliary table. */
2375 /*     UHERMT: Table of 2*(IORDRU+1) coeff. of 2*(IORDRU+1) polynoms of Hermite. */
2376 /*     VHERMT: Table of 2*(IORDRV+1) coeff. of 2*(IORDRV+1) polynoms of Hermite. */
2377 /*   SOSOTB: Table where the terms of constraints are added */
2378 /*           C(ui,vj) + C(ui,-vj) + C(-ui,vj) + C(-ui,-vj) */
2379 /*           with ui and vj positive roots of the Legendre polynom */
2380 /*           of degree NBPNTU and NBPNTV respectively. */
2381 /*   DISOTB: 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 polynom of Legendre */
2384 /*           of degree NBPNTU and NBPNTV respectively. */
2385 /*   SODITB: 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 /*   DIDITB: 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
2394 /*     COMMONS USED   : */
2395 /*     ---------------- */
2396
2397 /*     REFERENCES CALLED   : */
2398 /*     ----------------------- */
2399
2400 /*     DESCRIPTION/NOTES/LIMITATIONS : */
2401 /*     ----------------------------------- */
2402
2403 /* > */
2404 /* ********************************************************************** 
2405 */
2406
2407 /*   Name of the routine */
2408
2409
2410     /* Parameter adjustments */
2411     --urootl;
2412     diditb_dim1 = *nbpntu / 2 + 1;
2413     diditb_dim2 = *nbpntv / 2 + 1;
2414     diditb_offset = diditb_dim1 * diditb_dim2;
2415     diditb -= diditb_offset;
2416     disotb_dim1 = *nbpntu / 2;
2417     disotb_dim2 = *nbpntv / 2;
2418     disotb_offset = disotb_dim1 * (disotb_dim2 + 1) + 1;
2419     disotb -= disotb_offset;
2420     soditb_dim1 = *nbpntu / 2;
2421     soditb_dim2 = *nbpntv / 2;
2422     soditb_offset = soditb_dim1 * (soditb_dim2 + 1) + 1;
2423     soditb -= soditb_offset;
2424     sosotb_dim1 = *nbpntu / 2 + 1;
2425     sosotb_dim2 = *nbpntv / 2 + 1;
2426     sosotb_offset = sosotb_dim1 * sosotb_dim2;
2427     sosotb -= sosotb_offset;
2428     --vrootl;
2429     uhermt_dim1 = (*iordru << 1) + 2;
2430     uhermt_offset = uhermt_dim1;
2431     uhermt -= uhermt_offset;
2432     fpntbu_dim1 = *nbpntu;
2433     fpntbu_offset = fpntbu_dim1 + 1;
2434     fpntbu -= fpntbu_offset;
2435     vhermt_dim1 = (*iordrv << 1) + 2;
2436     vhermt_offset = vhermt_dim1;
2437     vhermt -= vhermt_offset;
2438     fpntbv_dim1 = *nbpntv;
2439     fpntbv_offset = fpntbv_dim1 + 1;
2440     fpntbv -= fpntbv_offset;
2441     contr4_dim1 = *ndimen;
2442     contr4_dim2 = *iordru + 2;
2443     contr4_offset = contr4_dim1 * (contr4_dim2 + 1) + 1;
2444     contr4 -= contr4_offset;
2445     contr3_dim1 = *ndimen;
2446     contr3_dim2 = *iordru + 2;
2447     contr3_offset = contr3_dim1 * (contr3_dim2 + 1) + 1;
2448     contr3 -= contr3_offset;
2449     contr2_dim1 = *ndimen;
2450     contr2_dim2 = *iordru + 2;
2451     contr2_offset = contr2_dim1 * (contr2_dim2 + 1) + 1;
2452     contr2 -= contr2_offset;
2453     contr1_dim1 = *ndimen;
2454     contr1_dim2 = *iordru + 2;
2455     contr1_offset = contr1_dim1 * (contr1_dim2 + 1) + 1;
2456     contr1 -= contr1_offset;
2457
2458     /* Function Body */
2459     ibb = AdvApp2Var_SysBase::mnfndeb_();
2460     if (ibb >= 3) {
2461         AdvApp2Var_SysBase::mgenmsg_("MMA2CD1", 7L);
2462     }
2463
2464 /* ------------------- Discretisation of Hermite polynoms ----------- 
2465 */
2466
2467     ncfhu = (*iordru + 1) << 1;
2468     i__1 = ncfhu;
2469     for (ii = 1; ii <= i__1; ++ii) {
2470         i__2 = *nbpntu;
2471         for (ll = 1; ll <= i__2; ++ll) {
2472             AdvApp2Var_MathBase::mmmpocur_(&ncfhu, &c__1, &ncfhu, &uhermt[ii * uhermt_dim1], &
2473                     urootl[ll], &fpntbu[ll + ii * fpntbu_dim1]);
2474 /* L20: */
2475         }
2476 /* L10: */
2477     }
2478     ncfhv = (*iordrv + 1) << 1;
2479     i__1 = ncfhv;
2480     for (jj = 1; jj <= i__1; ++jj) {
2481         i__2 = *nbpntv;
2482         for (kk = 1; kk <= i__2; ++kk) {
2483             AdvApp2Var_MathBase::mmmpocur_(&ncfhv, &c__1, &ncfhv, &vhermt[jj * vhermt_dim1], &
2484                     vrootl[kk], &fpntbv[kk + jj * fpntbv_dim1]);
2485 /* L40: */
2486         }
2487 /* L30: */
2488     }
2489
2490 /* ---- The discretizations of polynoms of constraints are subtracted ---- 
2491 */
2492
2493     nuroo = *nbpntu / 2;
2494     nvroo = *nbpntv / 2;
2495     i__1 = *ndimen;
2496     for (nd = 1; nd <= i__1; ++nd) {
2497
2498         i__2 = *iordrv + 1;
2499         for (jj = 1; jj <= i__2; ++jj) {
2500             i__3 = *iordru + 1;
2501             for (ii = 1; ii <= i__3; ++ii) {
2502                 bid1 = contr1[nd + (ii + jj * contr1_dim2) * contr1_dim1];
2503                 bid2 = contr2[nd + (ii + jj * contr2_dim2) * contr2_dim1];
2504                 bid3 = contr3[nd + (ii + jj * contr3_dim2) * contr3_dim1];
2505                 bid4 = contr4[nd + (ii + jj * contr4_dim2) * contr4_dim1];
2506
2507                 i__4 = nvroo;
2508                 for (kk = 1; kk <= i__4; ++kk) {
2509                     kkp = (*nbpntv + 1) / 2 + kk;
2510                     kkm = nvroo - kk + 1;
2511                     sov1 = fpntbv[kkp + ((jj << 1) - 1) * fpntbv_dim1] + 
2512                             fpntbv[kkm + ((jj << 1) - 1) * fpntbv_dim1];
2513                     div1 = fpntbv[kkp + ((jj << 1) - 1) * fpntbv_dim1] - 
2514                             fpntbv[kkm + ((jj << 1) - 1) * fpntbv_dim1];
2515                     sov2 = fpntbv[kkp + (jj << 1) * fpntbv_dim1] + fpntbv[kkm 
2516                             + (jj << 1) * fpntbv_dim1];
2517                     div2 = fpntbv[kkp + (jj << 1) * fpntbv_dim1] - fpntbv[kkm 
2518                             + (jj << 1) * fpntbv_dim1];
2519                     i__5 = nuroo;
2520                     for (ll = 1; ll <= i__5; ++ll) {
2521                         llp = (*nbpntu + 1) / 2 + ll;
2522                         llm = nuroo - ll + 1;
2523                         sou1 = fpntbu[llp + ((ii << 1) - 1) * fpntbu_dim1] + 
2524                                 fpntbu[llm + ((ii << 1) - 1) * fpntbu_dim1];
2525                         diu1 = fpntbu[llp + ((ii << 1) - 1) * fpntbu_dim1] - 
2526                                 fpntbu[llm + ((ii << 1) - 1) * fpntbu_dim1];
2527                         sou2 = fpntbu[llp + (ii << 1) * fpntbu_dim1] + fpntbu[
2528                                 llm + (ii << 1) * fpntbu_dim1];
2529                         diu2 = fpntbu[llp + (ii << 1) * fpntbu_dim1] - fpntbu[
2530                                 llm + (ii << 1) * fpntbu_dim1];
2531                         sosotb[ll + (kk + nd * sosotb_dim2) * sosotb_dim1] = 
2532                                 sosotb[ll + (kk + nd * sosotb_dim2) * 
2533                                 sosotb_dim1] - bid1 * sou1 * sov1 - bid2 * 
2534                                 sou2 * sov1 - bid3 * sou1 * sov2 - bid4 * 
2535                                 sou2 * sov2;
2536                         soditb[ll + (kk + nd * soditb_dim2) * soditb_dim1] = 
2537                                 soditb[ll + (kk + nd * soditb_dim2) * 
2538                                 soditb_dim1] - bid1 * sou1 * div1 - bid2 * 
2539                                 sou2 * div1 - bid3 * sou1 * div2 - bid4 * 
2540                                 sou2 * div2;
2541                         disotb[ll + (kk + nd * disotb_dim2) * disotb_dim1] = 
2542                                 disotb[ll + (kk + nd * disotb_dim2) * 
2543                                 disotb_dim1] - bid1 * diu1 * sov1 - bid2 * 
2544                                 diu2 * sov1 - bid3 * diu1 * sov2 - bid4 * 
2545                                 diu2 * sov2;
2546                         diditb[ll + (kk + nd * diditb_dim2) * diditb_dim1] = 
2547                                 diditb[ll + (kk + nd * diditb_dim2) * 
2548                                 diditb_dim1] - bid1 * diu1 * div1 - bid2 * 
2549                                 diu2 * div1 - bid3 * diu1 * div2 - bid4 * 
2550                                 diu2 * div2;
2551 /* L450: */
2552                     }
2553 /* L400: */
2554                 }
2555
2556 /* ------------ Case when the discretization is done only on the roots  
2557 ----------- */
2558 /* ----------   of Legendre polynom of uneven degree, 0 is root 
2559 ----------- */
2560
2561                 if (*nbpntu % 2 == 1) {
2562                     sou1 = fpntbu[nuroo + 1 + ((ii << 1) - 1) * fpntbu_dim1];
2563                     sou2 = fpntbu[nuroo + 1 + (ii << 1) * fpntbu_dim1];
2564                     i__4 = nvroo;
2565                     for (kk = 1; kk <= i__4; ++kk) {
2566                         kkp = (*nbpntv + 1) / 2 + kk;
2567                         kkm = nvroo - kk + 1;
2568                         sov1 = fpntbv[kkp + ((jj << 1) - 1) * fpntbv_dim1] + 
2569                                 fpntbv[kkm + ((jj << 1) - 1) * fpntbv_dim1];
2570                         div1 = fpntbv[kkp + ((jj << 1) - 1) * fpntbv_dim1] - 
2571                                 fpntbv[kkm + ((jj << 1) - 1) * fpntbv_dim1];
2572                         sov2 = fpntbv[kkp + (jj << 1) * fpntbv_dim1] + fpntbv[
2573                                 kkm + (jj << 1) * fpntbv_dim1];
2574                         div2 = fpntbv[kkp + (jj << 1) * fpntbv_dim1] - fpntbv[
2575                                 kkm + (jj << 1) * fpntbv_dim1];
2576                         sosotb[(kk + nd * sosotb_dim2) * sosotb_dim1] = 
2577                                 sosotb[(kk + nd * sosotb_dim2) * sosotb_dim1] 
2578                                 - bid1 * sou1 * sov1 - bid2 * sou2 * sov1 - 
2579                                 bid3 * sou1 * sov2 - bid4 * sou2 * sov2;
2580                         diditb[(kk + nd * diditb_dim2) * diditb_dim1] = 
2581                                 diditb[(kk + nd * diditb_dim2) * diditb_dim1] 
2582                                 - bid1 * sou1 * div1 - bid2 * sou2 * div1 - 
2583                                 bid3 * sou1 * div2 - bid4 * sou2 * div2;
2584 /* L500: */
2585                     }
2586                 }
2587
2588                 if (*nbpntv % 2 == 1) {
2589                     sov1 = fpntbv[nvroo + 1 + ((jj << 1) - 1) * fpntbv_dim1];
2590                     sov2 = fpntbv[nvroo + 1 + (jj << 1) * fpntbv_dim1];
2591                     i__4 = nuroo;
2592                     for (ll = 1; ll <= i__4; ++ll) {
2593                         llp = (*nbpntu + 1) / 2 + ll;
2594                         llm = nuroo - ll + 1;
2595                         sou1 = fpntbu[llp + ((ii << 1) - 1) * fpntbu_dim1] + 
2596                                 fpntbu[llm + ((ii << 1) - 1) * fpntbu_dim1];
2597                         diu1 = fpntbu[llp + ((ii << 1) - 1) * fpntbu_dim1] - 
2598                                 fpntbu[llm + ((ii << 1) - 1) * fpntbu_dim1];
2599                         sou2 = fpntbu[llp + (ii << 1) * fpntbu_dim1] + fpntbu[
2600                                 llm + (ii << 1) * fpntbu_dim1];
2601                         diu2 = fpntbu[llp + (ii << 1) * fpntbu_dim1] - fpntbu[
2602                                 llm + (ii << 1) * fpntbu_dim1];
2603                         sosotb[ll + nd * sosotb_dim2 * sosotb_dim1] = sosotb[
2604                                 ll + nd * sosotb_dim2 * sosotb_dim1] - bid1 * 
2605                                 sou1 * sov1 - bid2 * sou2 * sov1 - bid3 * 
2606                                 sou1 * sov2 - bid4 * sou2 * sov2;
2607                         diditb[ll + nd * diditb_dim2 * diditb_dim1] = diditb[
2608                                 ll + nd * diditb_dim2 * diditb_dim1] - bid1 * 
2609                                 diu1 * sov1 - bid2 * diu2 * sov1 - bid3 * 
2610                                 diu1 * sov2 - bid4 * diu2 * sov2;
2611 /* L600: */
2612                     }
2613                 }
2614
2615                 if (*nbpntu % 2 == 1 && *nbpntv % 2 == 1) {
2616                     sou1 = fpntbu[nuroo + 1 + ((ii << 1) - 1) * fpntbu_dim1];
2617                     sou2 = fpntbu[nuroo + 1 + (ii << 1) * fpntbu_dim1];
2618                     sov1 = fpntbv[nvroo + 1 + ((jj << 1) - 1) * fpntbv_dim1];
2619                     sov2 = fpntbv[nvroo + 1 + (jj << 1) * fpntbv_dim1];
2620                     sosotb[nd * sosotb_dim2 * sosotb_dim1] = sosotb[nd * 
2621                             sosotb_dim2 * sosotb_dim1] - bid1 * sou1 * sov1 - 
2622                             bid2 * sou2 * sov1 - bid3 * sou1 * sov2 - bid4 * 
2623                             sou2 * sov2;
2624                     diditb[nd * diditb_dim2 * diditb_dim1] = diditb[nd * 
2625                             diditb_dim2 * diditb_dim1] - bid1 * sou1 * sov1 - 
2626                             bid2 * sou2 * sov1 - bid3 * sou1 * sov2 - bid4 * 
2627                             sou2 * sov2;
2628                 }
2629
2630 /* L300: */
2631             }
2632 /* L200: */
2633         }
2634 /* L100: */
2635     }
2636     goto L9999;
2637
2638 /* ------------------------------ The End ------------------------------- 
2639 */
2640
2641 L9999:
2642     if (ibb >= 3) {
2643         AdvApp2Var_SysBase::mgsomsg_("MMA2CD1", 7L);
2644     }
2645     return 0;
2646 } /* mma2cd1_ */
2647
2648 //=======================================================================
2649 //function : mma2cd2_
2650 //purpose  : 
2651 //=======================================================================
2652 int mma2cd2_(integer *ndimen, 
2653              integer *nbpntu, 
2654              integer *nbpntv, 
2655              doublereal *vrootl, 
2656              integer *iordrv, 
2657              doublereal *sotbv1, 
2658              doublereal *sotbv2, 
2659              doublereal *ditbv1, 
2660              doublereal *ditbv2, 
2661              doublereal *fpntab, 
2662              doublereal *vhermt, 
2663              doublereal *sosotb, 
2664              doublereal *soditb, 
2665              doublereal *disotb, 
2666              doublereal *diditb)
2667
2668 {
2669   static integer c__1 = 1;
2670   /* System generated locals */
2671   integer sotbv1_dim1, sotbv1_dim2, sotbv1_offset, sotbv2_dim1, sotbv2_dim2,
2672   sotbv2_offset, ditbv1_dim1, ditbv1_dim2, ditbv1_offset, 
2673   ditbv2_dim1, ditbv2_dim2, ditbv2_offset, fpntab_dim1, 
2674   fpntab_offset, vhermt_dim1, vhermt_offset, sosotb_dim1, 
2675   sosotb_dim2, sosotb_offset, diditb_dim1, diditb_dim2, 
2676   diditb_offset, soditb_dim1, soditb_dim2, soditb_offset, 
2677   disotb_dim1, disotb_dim2, disotb_offset, i__1, i__2, i__3, i__4;
2678
2679   /* Local variables */
2680   static integer ncfhv, nuroo, nvroo, ii, nd, jj, kk, ibb, jjm, jjp;
2681   static doublereal bid1, bid2, bid3, bid4;
2682
2683 /* ********************************************************************** 
2684 */
2685 /*     FUNCTION : */
2686 /*     ---------- */
2687 /*     Discretisation on the parameters of polynoms of interpolation */
2688 /*     of constraints on 2 borders iso-V of order IORDRV. */
2689
2690
2691 /*     KEYWORDS : */
2692 /*     ----------- */
2693 /*     TOUS, AB_SPECIFI::CONTRAINTE&, DISCRETISATION, &POINT */
2694
2695
2696
2697 /*     INPUT ARGUMENTS : */
2698 /*     ------------------ */
2699 /*     NDIMEN: Dimension of the space. */
2700 /*     NBPNTU: Nb of INTERNAL parameters of discretisation by U. */
2701 /*             This is also the nb of root of Legendre polynom where discretization is done. */
2702 /*     UROOTL: Table of parameters of discretisation ON (-1,1) by U. 
2703 */
2704 /*     NBPNTV: Nb of INTERNAL  parameters of discretisation by V. */
2705 /*             This is also the nb of root of Legendre polynom where discretization is done. */
2706 /*     VROOTL: Table of discretization parameters on (-1,1) by V. 
2707 /*     IORDRV: Order of constraint imposed at the extremities of iso-V */
2708 /*             = 0, calculate the extremities of iso-V */
2709 /*             = 1, calculate, additionally, the 1st derivative in the direction of iso-V */
2710 /*             = 2, calculate, additionally, the 2nd derivative in the direction of iso-V */
2711 /*     SOTBV1: Table of NBPNTV/2 sums of 2 index points  */
2712 /*             NBPNTV-II+1 and II, for II = 1, NBPNTV/2 on iso-V0. */
2713 /*     SOTBV2: Table of NBPNTV/2 sums of 2 index points  */
2714 /*             NBPNTV-II+1 and II, for II = 1, NBPNTV/2 on iso-V1. */
2715 /*     DITBV1: Table of NBPNTV/2 differences of 2 index points */
2716 /*             NBPNTV-II+1 and II, for II = 1, NBPNTV/2 on iso-V0. */
2717 /*     DITBV2: Table of NBPNTV/2 differences of 2 index points */
2718 /*             NBPNTV-II+1 and II, for II = 1, NBPNTV/2 on iso-V1. */
2719 /*     SOSOTB: Preinitialized table (input/output argument). */
2720 /*     DISOTB: Preinitialized table (input/output argument). */
2721 /*     SODITB: Preinitialized table (input/output argument). */
2722 /*     DIDITB: Preinitialized table (input/output argument) */
2723
2724 /*     OUTPUT ARGUMENTS : */
2725 /*     ------------------- */
2726 /*     FPNTAB: Auxiliary table. */
2727 /*     VHERMT: Table of 2*(IORDRV+1) coeff. of 2*(IORDRV+1) polynoms of Hermite. */
2728 /*   SOSOTB: Table where the terms of constraints are added */
2729 /*           C(ui,vj) + C(ui,-vj) + C(-ui,vj) + C(-ui,-vj) */
2730 /*           with ui and vj positive roots of the Legendre polynom */
2731 /*           of degree NBPNTU and NBPNTV respectively. */
2732 /*   DISOTB: 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 polynom of Legendre */
2735 /*           of degree NBPNTU and NBPNTV respectively. */
2736 /*   SODITB: 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 /*   DIDITB: 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
2745 /*     COMMONS USED   : */
2746 /*     ---------------- */
2747
2748 /*     REFERENCES CALLED   : */
2749 /*     ----------------------- */
2750
2751 /*     DESCRIPTION/NOTES/LIMITATIONS : */
2752 /*     ----------------------------------- */
2753
2754
2755 /* > */
2756 /* ********************************************************************** 
2757 */
2758
2759 /*   Name of the routine */
2760
2761
2762     /* Parameter adjustments */
2763     diditb_dim1 = *nbpntu / 2 + 1;
2764     diditb_dim2 = *nbpntv / 2 + 1;
2765     diditb_offset = diditb_dim1 * diditb_dim2;
2766     diditb -= diditb_offset;
2767     disotb_dim1 = *nbpntu / 2;
2768     disotb_dim2 = *nbpntv / 2;
2769     disotb_offset = disotb_dim1 * (disotb_dim2 + 1) + 1;
2770     disotb -= disotb_offset;
2771     soditb_dim1 = *nbpntu / 2;
2772     soditb_dim2 = *nbpntv / 2;
2773     soditb_offset = soditb_dim1 * (soditb_dim2 + 1) + 1;
2774     soditb -= soditb_offset;
2775     sosotb_dim1 = *nbpntu / 2 + 1;
2776     sosotb_dim2 = *nbpntv / 2 + 1;
2777     sosotb_offset = sosotb_dim1 * sosotb_dim2;
2778     sosotb -= sosotb_offset;
2779     --vrootl;
2780     vhermt_dim1 = (*iordrv << 1) + 2;
2781     vhermt_offset = vhermt_dim1;
2782     vhermt -= vhermt_offset;
2783     fpntab_dim1 = *nbpntv;
2784     fpntab_offset = fpntab_dim1 + 1;
2785     fpntab -= fpntab_offset;
2786     ditbv2_dim1 = *nbpntu / 2 + 1;
2787     ditbv2_dim2 = *ndimen;
2788     ditbv2_offset = ditbv2_dim1 * (ditbv2_dim2 + 1);
2789     ditbv2 -= ditbv2_offset;
2790     ditbv1_dim1 = *nbpntu / 2 + 1;
2791     ditbv1_dim2 = *ndimen;
2792     ditbv1_offset = ditbv1_dim1 * (ditbv1_dim2 + 1);
2793     ditbv1 -= ditbv1_offset;
2794     sotbv2_dim1 = *nbpntu / 2 + 1;
2795     sotbv2_dim2 = *ndimen;
2796     sotbv2_offset = sotbv2_dim1 * (sotbv2_dim2 + 1);
2797     sotbv2 -= sotbv2_offset;
2798     sotbv1_dim1 = *nbpntu / 2 + 1;
2799     sotbv1_dim2 = *ndimen;
2800     sotbv1_offset = sotbv1_dim1 * (sotbv1_dim2 + 1);
2801     sotbv1 -= sotbv1_offset;
2802
2803     /* Function Body */
2804     ibb = AdvApp2Var_SysBase::mnfndeb_();
2805     if (ibb >= 3) {
2806         AdvApp2Var_SysBase::mgenmsg_("MMA2CD2", 7L);
2807     }
2808
2809 /* ------------------- Discretization of Hermit polynoms ----------- 
2810 */
2811
2812     ncfhv = (*iordrv + 1) << 1;
2813     i__1 = ncfhv;
2814     for (ii = 1; ii <= i__1; ++ii) {
2815         i__2 = *nbpntv;
2816         for (jj = 1; jj <= i__2; ++jj) {
2817             AdvApp2Var_MathBase::mmmpocur_(&ncfhv, &c__1, &ncfhv, &vhermt[ii * vhermt_dim1], &
2818                     vrootl[jj], &fpntab[jj + ii * fpntab_dim1]);
2819 /* L60: */
2820         }
2821 /* L50: */
2822     }
2823
2824 /* ---- The discretizations of polynoms of constraints are subtracted ---- 
2825 */
2826
2827     nuroo = *nbpntu / 2;
2828     nvroo = *nbpntv / 2;
2829
2830     i__1 = *ndimen;
2831     for (nd = 1; nd <= i__1; ++nd) {
2832         i__2 = *iordrv + 1;
2833         for (ii = 1; ii <= i__2; ++ii) {
2834
2835             i__3 = nuroo;
2836             for (kk = 1; kk <= i__3; ++kk) {
2837                 bid1 = sotbv1[kk + (nd + ii * sotbv1_dim2) * sotbv1_dim1];
2838                 bid2 = sotbv2[kk + (nd + ii * sotbv2_dim2) * sotbv2_dim1];
2839                 bid3 = ditbv1[kk + (nd + ii * ditbv1_dim2) * ditbv1_dim1];
2840                 bid4 = ditbv2[kk + (nd + ii * ditbv2_dim2) * ditbv2_dim1];
2841                 i__4 = nvroo;
2842                 for (jj = 1; jj <= i__4; ++jj) {
2843                     jjp = (*nbpntv + 1) / 2 + jj;
2844                     jjm = nvroo - jj + 1;
2845                     sosotb[kk + (jj + nd * sosotb_dim2) * sosotb_dim1] = 
2846                             sosotb[kk + (jj + nd * sosotb_dim2) * sosotb_dim1]
2847                              - bid1 * (fpntab[jjp + ((ii << 1) - 1) * 
2848                             fpntab_dim1] + fpntab[jjm + ((ii << 1) - 1) * 
2849                             fpntab_dim1]) - bid2 * (fpntab[jjp + (ii << 1) * 
2850                             fpntab_dim1] + fpntab[jjm + (ii << 1) * 
2851                             fpntab_dim1]);
2852                     disotb[kk + (jj + nd * disotb_dim2) * disotb_dim1] = 
2853                             disotb[kk + (jj + nd * disotb_dim2) * disotb_dim1]
2854                              - bid3 * (fpntab[jjp + ((ii << 1) - 1) * 
2855                             fpntab_dim1] + fpntab[jjm + ((ii << 1) - 1) * 
2856                             fpntab_dim1]) - bid4 * (fpntab[jjp + (ii << 1) * 
2857                             fpntab_dim1] + fpntab[jjm + (ii << 1) * 
2858                             fpntab_dim1]);
2859                     soditb[kk + (jj + nd * soditb_dim2) * soditb_dim1] = 
2860                             soditb[kk + (jj + nd * soditb_dim2) * soditb_dim1]
2861                              - bid1 * (fpntab[jjp + ((ii << 1) - 1) * 
2862                             fpntab_dim1] - fpntab[jjm + ((ii << 1) - 1) * 
2863                             fpntab_dim1]) - bid2 * (fpntab[jjp + (ii << 1) * 
2864                             fpntab_dim1] - fpntab[jjm + (ii << 1) * 
2865                             fpntab_dim1]);
2866                     diditb[kk + (jj + nd * diditb_dim2) * diditb_dim1] = 
2867                             diditb[kk + (jj + nd * diditb_dim2) * diditb_dim1]
2868                              - bid3 * (fpntab[jjp + ((ii << 1) - 1) * 
2869                             fpntab_dim1] - fpntab[jjm + ((ii << 1) - 1) * 
2870                             fpntab_dim1]) - bid4 * (fpntab[jjp + (ii << 1) * 
2871                             fpntab_dim1] - fpntab[jjm + (ii << 1) * 
2872                             fpntab_dim1]);
2873 /* L400: */
2874                 }
2875 /* L300: */
2876             }
2877 /* L200: */
2878         }
2879
2880 /* ------------ Case when the discretization is done only on the roots  */
2881 /* ----------   of Legendre polynom of uneven degree, 0 is root */
2882
2883
2884         if (*nbpntv % 2 == 1) {
2885             i__2 = *iordrv + 1;
2886             for (ii = 1; ii <= i__2; ++ii) {
2887                 i__3 = nuroo;
2888                 for (kk = 1; kk <= i__3; ++kk) {
2889                     bid1 = sotbv1[kk + (nd + ii * sotbv1_dim2) * sotbv1_dim1] 
2890                             * fpntab[nvroo + 1 + ((ii << 1) - 1) * 
2891                             fpntab_dim1] + sotbv2[kk + (nd + ii * sotbv2_dim2)
2892                              * sotbv2_dim1] * fpntab[nvroo + 1 + (ii << 1) * 
2893                             fpntab_dim1];
2894                     sosotb[kk + nd * sosotb_dim2 * sosotb_dim1] -= bid1;
2895                     bid2 = ditbv1[kk + (nd + ii * ditbv1_dim2) * ditbv1_dim1] 
2896                             * fpntab[nvroo + 1 + ((ii << 1) - 1) * 
2897                             fpntab_dim1] + ditbv2[kk + (nd + ii * ditbv2_dim2)
2898                              * ditbv2_dim1] * fpntab[nvroo + 1 + (ii << 1) * 
2899                             fpntab_dim1];
2900                     diditb[kk + nd * diditb_dim2 * diditb_dim1] -= bid2;
2901 /* L550: */
2902                 }
2903 /* L500: */
2904             }
2905         }
2906
2907         if (*nbpntu % 2 == 1) {
2908             i__2 = *iordrv + 1;
2909             for (ii = 1; ii <= i__2; ++ii) {
2910                 i__3 = nvroo;
2911                 for (jj = 1; jj <= i__3; ++jj) {
2912                     jjp = (*nbpntv + 1) / 2 + jj;
2913                     jjm = nvroo - jj + 1;
2914                     bid1 = sotbv1[(nd + ii * sotbv1_dim2) * sotbv1_dim1] * (
2915                             fpntab[jjp + ((ii << 1) - 1) * fpntab_dim1] + 
2916                             fpntab[jjm + ((ii << 1) - 1) * fpntab_dim1]) + 
2917                             sotbv2[(nd + ii * sotbv2_dim2) * sotbv2_dim1] * (
2918                             fpntab[jjp + (ii << 1) * fpntab_dim1] + fpntab[
2919                             jjm + (ii << 1) * fpntab_dim1]);
2920                     sosotb[(jj + nd * sosotb_dim2) * sosotb_dim1] -= bid1;
2921                     bid2 = sotbv1[(nd + ii * sotbv1_dim2) * sotbv1_dim1] * (
2922                             fpntab[jjp + ((ii << 1) - 1) * fpntab_dim1] - 
2923                             fpntab[jjm + ((ii << 1) - 1) * fpntab_dim1]) + 
2924                             sotbv2[(nd + ii * sotbv2_dim2) * sotbv2_dim1] * (
2925                             fpntab[jjp + (ii << 1) * fpntab_dim1] - fpntab[
2926                             jjm + (ii << 1) * fpntab_dim1]);
2927                     diditb[jj + nd * diditb_dim2 * diditb_dim1] -= bid2;
2928 /* L650: */
2929                 }
2930 /* L600: */
2931             }
2932         }
2933
2934         if (*nbpntu % 2 == 1 && *nbpntv % 2 == 1) {
2935             i__2 = *iordrv + 1;
2936             for (ii = 1; ii <= i__2; ++ii) {
2937                 bid1 = sotbv1[(nd + ii * sotbv1_dim2) * sotbv1_dim1] * fpntab[
2938                         nvroo + 1 + ((ii << 1) - 1) * fpntab_dim1] + sotbv2[(
2939                         nd + ii * sotbv2_dim2) * sotbv2_dim1] * fpntab[nvroo 
2940                         + 1 + (ii << 1) * fpntab_dim1];
2941                 sosotb[nd * sosotb_dim2 * sosotb_dim1] -= bid1;
2942 /* L700: */
2943             }
2944         }
2945
2946 /* L100: */
2947     }
2948     goto L9999;
2949
2950 /* ------------------------------ The End ------------------------------- 
2951 */
2952
2953 L9999:
2954     if (ibb >= 3) {
2955         AdvApp2Var_SysBase::mgsomsg_("MMA2CD2", 7L);
2956     }
2957     return 0;
2958 } /* mma2cd2_ */
2959
2960 //=======================================================================
2961 //function : mma2cd3_
2962 //purpose  : 
2963 //=======================================================================
2964 int mma2cd3_(integer *ndimen,
2965              integer *nbpntu, 
2966              doublereal *urootl, 
2967              integer *nbpntv, 
2968              integer *iordru, 
2969              doublereal *sotbu1, 
2970              doublereal *sotbu2, 
2971              doublereal *ditbu1, 
2972              doublereal *ditbu2, 
2973              doublereal *fpntab, 
2974              doublereal *uhermt, 
2975              doublereal *sosotb, 
2976              doublereal *soditb, 
2977              doublereal *disotb, 
2978              doublereal *diditb)
2979
2980 {
2981   static integer c__1 = 1;
2982
2983    /* System generated locals */
2984     integer sotbu1_dim1, sotbu1_dim2, sotbu1_offset, sotbu2_dim1, sotbu2_dim2,
2985              sotbu2_offset, ditbu1_dim1, ditbu1_dim2, ditbu1_offset, 
2986             ditbu2_dim1, ditbu2_dim2, ditbu2_offset, fpntab_dim1, 
2987             fpntab_offset, uhermt_dim1, uhermt_offset, sosotb_dim1, 
2988             sosotb_dim2, sosotb_offset, diditb_dim1, diditb_dim2, 
2989             diditb_offset, soditb_dim1, soditb_dim2, soditb_offset, 
2990             disotb_dim1, disotb_dim2, disotb_offset, i__1, i__2, i__3, i__4;
2991
2992     /* Local variables */
2993     static integer ncfhu, nuroo, nvroo, ii, nd, jj, kk, ibb, kkm, kkp;
2994     static doublereal bid1, bid2, bid3, bid4;
2995
2996 /* ********************************************************************** 
2997 */
2998 /*     FUNCTION : */
2999 /*     ---------- */
3000 /*     Discretisation on the parameters of polynoms of interpolation */
3001 /*     of constraints on 2 borders iso-U of order IORDRU. */
3002
3003
3004 /*     KEYWORDS : */
3005 /*     ----------- */
3006 /*     TOUS, AB_SPECIFI::CONTRAINTE&, DISCRETISATION, &POINT */
3007
3008 /*     INPUT ARGUMENTS : */
3009 /*     ------------------ */
3010 /*     NDIMEN: Dimension of the space. */
3011 /*     NBPNTU: Nb of INTERNAL parameters of discretisation by U. */
3012 /*             This is also the nb of root of Legendre polynom where discretization is done. */
3013 /*     UROOTL: Table of parameters of discretisation ON (-1,1) by U. 
3014 */
3015 /*     NBPNTV: Nb of INTERNAL  parameters of discretisation by V. */
3016 /*             This is also the nb of root of Legendre polynom where discretization is done. */
3017 /*     IORDRV: Order of constraint imposed at the extremities of iso-V */
3018 /*             = 0, calculate the extremities of iso-V */
3019 /*             = 1, calculate, additionally, the 1st derivative in the direction of iso-V */
3020 /*             = 2, calculate, additionally, the 2nd derivative in the direction of iso-V */
3021 /*     SOTBU1: Table of NBPNTU/2 sums of 2 index points  */
3022 /*             NBPNTU-II+1 and II, for II = 1, NBPNTU/2 on iso-V0. */
3023 /*     SOTBU2: Table of NBPNTV/2 sums of 2 index points  */
3024 /*             NBPNTU-II+1 and II, for II = 1, NBPNTU/2 on iso-V1. */
3025 /*     DITBU1: Table of NBPNTU/2 differences of 2 index points */
3026 /*             NBPNTU-II+1 and II, for II = 1, NBPNTU/2 on iso-V0. */
3027 /*     DITBU2: Table of NBPNTU/2 differences of 2 index points */
3028 /*             NBPNTU-II+1 and II, for II = 1, NBPNTU/2 on iso-V1. */
3029 /*     SOSOTB: Preinitialized table (input/output argument). */
3030 /*     DISOTB: Preinitialized table (input/output argument). */
3031 /*     SODITB: Preinitialized table (input/output argument). */
3032 /*     DIDITB: Preinitialized table (input/output argument) */
3033
3034 /*     OUTPUT ARGUMENTS : */
3035 /*     ------------------- */
3036 /*     FPNTAB: Auxiliary table. */
3037 /*     UHERMT: Table of 2*(IORDRU+1) coeff. of 2*(IORDRU+1) polynoms of Hermite. */
3038 /*   SOSOTB: Table where the terms of constraints are added */
3039 /*           C(ui,vj) + C(ui,-vj) + C(-ui,vj) + C(-ui,-vj) */
3040 /*           with ui and vj positive roots of the Legendre polynom */
3041 /*           of degree NBPNTU and NBPNTV respectively. */
3042 /*   DISOTB: 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 polynom of Legendre */
3045 /*           of degree NBPNTU and NBPNTV respectively. */
3046 /*   SODITB: 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 /*   DIDITB: 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
3055 /*     COMMONS USED   : */
3056 /*     ---------------- */
3057
3058 /*     REFERENCES CALLED   : */
3059 /*     ----------------------- */
3060
3061 /*     DESCRIPTION/NOTES/LIMITATIONS : */
3062 /*     ----------------------------------- */
3063
3064 /* $    HISTORIQUE DES MODIFICATIONS   : */
3065 /*     -------------------------------- */
3066 /*     08-08-1991: RBD; Creation. */
3067 /* > */
3068 /* ********************************************************************** 
3069 */
3070
3071 /*   Name of the routine */
3072
3073
3074     /* Parameter adjustments */
3075     --urootl;
3076     diditb_dim1 = *nbpntu / 2 + 1;
3077     diditb_dim2 = *nbpntv / 2 + 1;
3078     diditb_offset = diditb_dim1 * diditb_dim2;
3079     diditb -= diditb_offset;
3080     disotb_dim1 = *nbpntu / 2;
3081     disotb_dim2 = *nbpntv / 2;
3082     disotb_offset = disotb_dim1 * (disotb_dim2 + 1) + 1;
3083     disotb -= disotb_offset;
3084     soditb_dim1 = *nbpntu / 2;
3085     soditb_dim2 = *nbpntv / 2;
3086     soditb_offset = soditb_dim1 * (soditb_dim2 + 1) + 1;
3087     soditb -= soditb_offset;
3088     sosotb_dim1 = *nbpntu / 2 + 1;
3089     sosotb_dim2 = *nbpntv / 2 + 1;
3090     sosotb_offset = sosotb_dim1 * sosotb_dim2;
3091     sosotb -= sosotb_offset;
3092     uhermt_dim1 = (*iordru << 1) + 2;
3093     uhermt_offset = uhermt_dim1;
3094     uhermt -= uhermt_offset;
3095     fpntab_dim1 = *nbpntu;
3096     fpntab_offset = fpntab_dim1 + 1;
3097     fpntab -= fpntab_offset;
3098     ditbu2_dim1 = *nbpntv / 2 + 1;
3099     ditbu2_dim2 = *ndimen;
3100     ditbu2_offset = ditbu2_dim1 * (ditbu2_dim2 + 1);
3101     ditbu2 -= ditbu2_offset;
3102     ditbu1_dim1 = *nbpntv / 2 + 1;
3103     ditbu1_dim2 = *ndimen;
3104     ditbu1_offset = ditbu1_dim1 * (ditbu1_dim2 + 1);
3105     ditbu1 -= ditbu1_offset;
3106     sotbu2_dim1 = *nbpntv / 2 + 1;
3107     sotbu2_dim2 = *ndimen;
3108     sotbu2_offset = sotbu2_dim1 * (sotbu2_dim2 + 1);
3109     sotbu2 -= sotbu2_offset;
3110     sotbu1_dim1 = *nbpntv / 2 + 1;
3111     sotbu1_dim2 = *ndimen;
3112     sotbu1_offset = sotbu1_dim1 * (sotbu1_dim2 + 1);
3113     sotbu1 -= sotbu1_offset;
3114
3115     /* Function Body */
3116     ibb = AdvApp2Var_SysBase::mnfndeb_();
3117     if (ibb >= 3) {
3118         AdvApp2Var_SysBase::mgenmsg_("MMA2CD3", 7L);
3119     }
3120
3121 /* ------------------- Discretization of polynoms of Hermit ----------- 
3122 */
3123
3124     ncfhu = (*iordru + 1) << 1;
3125     i__1 = ncfhu;
3126     for (ii = 1; ii <= i__1; ++ii) {
3127         i__2 = *nbpntu;
3128         for (kk = 1; kk <= i__2; ++kk) {
3129             AdvApp2Var_MathBase::mmmpocur_(&ncfhu, 
3130                                            &c__1, 
3131                                            &ncfhu,
3132                                            &uhermt[ii * uhermt_dim1],
3133                                            &urootl[kk], 
3134                                            &fpntab[kk + ii * fpntab_dim1]);
3135 /* L60: */
3136         }
3137 /* L50: */
3138     }
3139
3140 /* ---- The discretizations of polynoms of constraints are subtracted ---- 
3141 */
3142
3143     nvroo = *nbpntv / 2;
3144     nuroo = *nbpntu / 2;
3145
3146     i__1 = *ndimen;
3147     for (nd = 1; nd <= i__1; ++nd) {
3148         i__2 = *iordru + 1;
3149         for (ii = 1; ii <= i__2; ++ii) {
3150
3151             i__3 = nvroo;
3152             for (jj = 1; jj <= i__3; ++jj) {
3153                 bid1 = sotbu1[jj + (nd + ii * sotbu1_dim2) * sotbu1_dim1];
3154                 bid2 = sotbu2[jj + (nd + ii * sotbu2_dim2) * sotbu2_dim1];
3155                 bid3 = ditbu1[jj + (nd + ii * ditbu1_dim2) * ditbu1_dim1];
3156                 bid4 = ditbu2[jj + (nd + ii * ditbu2_dim2) * ditbu2_dim1];
3157                 i__4 = nuroo;
3158                 for (kk = 1; kk <= i__4; ++kk) {
3159                     kkp = (*nbpntu + 1) / 2 + kk;
3160                     kkm = nuroo - kk + 1;
3161                     sosotb[kk + (jj + nd * sosotb_dim2) * sosotb_dim1] = 
3162                             sosotb[kk + (jj + nd * sosotb_dim2) * sosotb_dim1]
3163                              - bid1 * (fpntab[kkp + ((ii << 1) - 1) * 
3164                             fpntab_dim1] + fpntab[kkm + ((ii << 1) - 1) * 
3165                             fpntab_dim1]) - bid2 * (fpntab[kkp + (ii << 1) * 
3166                             fpntab_dim1] + fpntab[kkm + (ii << 1) * 
3167                             fpntab_dim1]);
3168                     disotb[kk + (jj + nd * disotb_dim2) * disotb_dim1] = 
3169                             disotb[kk + (jj + nd * disotb_dim2) * disotb_dim1]
3170                              - bid1 * (fpntab[kkp + ((ii << 1) - 1) * 
3171                             fpntab_dim1] - fpntab[kkm + ((ii << 1) - 1) * 
3172                             fpntab_dim1]) - bid2 * (fpntab[kkp + (ii << 1) * 
3173                             fpntab_dim1] - fpntab[kkm + (ii << 1) * 
3174                             fpntab_dim1]);
3175                     soditb[kk + (jj + nd * soditb_dim2) * soditb_dim1] = 
3176                             soditb[kk + (jj + nd * soditb_dim2) * soditb_dim1]
3177                              - bid3 * (fpntab[kkp + ((ii << 1) - 1) * 
3178                             fpntab_dim1] + fpntab[kkm + ((ii << 1) - 1) * 
3179                             fpntab_dim1]) - bid4 * (fpntab[kkp + (ii << 1) * 
3180                             fpntab_dim1] + fpntab[kkm + (ii << 1) * 
3181                             fpntab_dim1]);
3182                     diditb[kk + (jj + nd * diditb_dim2) * diditb_dim1] = 
3183                             diditb[kk + (jj + nd * diditb_dim2) * diditb_dim1]
3184                              - bid3 * (fpntab[kkp + ((ii << 1) - 1) * 
3185                             fpntab_dim1] - fpntab[kkm + ((ii << 1) - 1) * 
3186                             fpntab_dim1]) - bid4 * (fpntab[kkp + (ii << 1) * 
3187                             fpntab_dim1] - fpntab[kkm + (ii << 1) * 
3188                             fpntab_dim1]);
3189 /* L400: */
3190                 }
3191 /* L300: */
3192             }
3193 /* L200: */
3194         }
3195
3196 /* ------------ Case when the discretization is done only on the roots  */
3197 /* ----------   of Legendre polynom of uneven degree, 0 is root */
3198
3199
3200
3201         if (*nbpntu % 2 == 1) {
3202             i__2 = *iordru + 1;
3203             for (ii = 1; ii <= i__2; ++ii) {
3204                 i__3 = nvroo;
3205                 for (jj = 1; jj <= i__3; ++jj) {
3206                     bid1 = sotbu1[jj + (nd + ii * sotbu1_dim2) * sotbu1_dim1] 
3207                             * fpntab[nuroo + 1 + ((ii << 1) - 1) * 
3208                             fpntab_dim1] + sotbu2[jj + (nd + ii * sotbu2_dim2)
3209                              * sotbu2_dim1] * fpntab[nuroo + 1 + (ii << 1) * 
3210                             fpntab_dim1];
3211                     sosotb[(jj + nd * sosotb_dim2) * sosotb_dim1] -= bid1;
3212                     bid2 = ditbu1[jj + (nd + ii * ditbu1_dim2) * ditbu1_dim1] 
3213                             * fpntab[nuroo + 1 + ((ii << 1) - 1) * 
3214                             fpntab_dim1] + ditbu2[jj + (nd + ii * ditbu2_dim2)
3215                              * ditbu2_dim1] * fpntab[nuroo + 1 + (ii << 1) * 
3216                             fpntab_dim1];
3217                     diditb[(jj + nd * diditb_dim2) * diditb_dim1] -= bid2;