2 // AdvApp2Var_ApproxF2var.cxx
5 #include <AdvApp2Var_SysBase.hxx>
6 #include <AdvApp2Var_MathBase.hxx>
7 #include <AdvApp2Var_Data_f2c.hxx>
8 #include <AdvApp2Var_Data.hxx>
9 #include <AdvApp2Var_ApproxF2var.hxx>
13 int mmjacpt_(const integer *ndimen,
14 const integer *ncoefu,
15 const integer *ncoefv,
16 const integer *iordru,
17 const integer *iordrv,
18 const doublereal *ptclgd,
25 int mma2ce2_(integer *numdec,
60 int mma2cfu_(integer *ndujac,
72 int mma2cfv_(integer *ndvjac,
82 int mma2er1_(integer *ndjacu,
98 int mma2er2_(integer *ndjacu,
117 int mma2moy_(integer *ndgumx,
130 int mma2ds2_(integer *ndimen,
133 const AdvApp2Var_EvaluatorFunc2Var& foncnp,
151 int mma1fdi_(integer *ndimen,
153 const AdvApp2Var_EvaluatorFunc2Var& foncnp,
168 int mma1cdi_(integer *ndimen,
180 int mma1jak_(integer *ndimen,
190 int mma1cnt_(integer *ndimen,
199 int mma1fer_(integer *ndimen,
214 int mma1noc_(doublereal *dfuvin,
225 int mmmapcoe_(integer *ndim,
235 int mmaperm_(integer *ncofmx,
244 #define mmapgss_1 mmapgss_
245 #define mmapgs0_1 mmapgs0_
246 #define mmapgs1_1 mmapgs1_
247 #define mmapgs2_1 mmapgs2_
249 //=======================================================================
250 //function : mma1cdi_
252 //=======================================================================
253 int mma1cdi_(integer *ndimen,
265 static integer c__1 = 1;
267 /* System generated locals */
268 integer contr1_dim1, contr1_offset, contr2_dim1, contr2_offset,
269 somtab_dim1, somtab_offset, diftab_dim1, diftab_offset,
270 fpntab_dim1, fpntab_offset, hermit_dim1, hermit_offset, i__1,
273 /* Local variables */
274 static integer nroo2, ncfhe, nd, ii, kk;
275 static integer ibb, kkm, kkp;
276 static doublereal bid1, bid2, bid3;
278 /* **********************************************************************
282 /* Discretisation on the parameters of interpolation polynomes */
283 /* constraints of order IORDRE. */
287 /* ALL, AB_SPECIFI::CONTRAINTE&, DISCRETISATION, &POINT */
289 /* INPUT ARGUMENTS : */
290 /* ------------------ */
291 /* NDIMEN: Space dimension. */
292 /* NBROOT: Number of INTERNAL discretisation parameters. */
293 /* It is also the root number Legendre polynome where */
294 /* the discretization is performed. */
295 /* ROOTLG: Table of discretization parameters ON (-1,1). */
296 /* IORDRE: Order of constraint imposed to the extremities of the iso. */
297 /* = 0, the extremities of the iso are calculated */
298 /* = 1, additionally, the 1st derivative in the direction */
299 /* of the iso is calculated. */
300 /* = 2, additionally, the 2nd derivative in the direction */
301 /* of the iso is calculated. */
302 /* CONTR1: Contains, if IORDRE>=0, values IORDRE+1 in TTABLE(0)
304 /* (1st extremity) of derivatives of F(Uc,Ve) or F(Ue,Vc), */
306 /* CONTR2: Contains, if IORDRE>=0, values IORDRE+1 in */
307 /* TTABLE(NBROOT+1) (2nd extremity) of: */
308 /* If ISOFAV=1, derived of order IDERIV by U, derived */
309 /* ordre 0 to IORDRE by V of F(Uc,Ve) or Uc=TCONST */
310 /* (fixed iso value) and Ve is the fixed extremity. */
311 /* If ISOFAV=2, derivative of order IDERIV by V, derivative */
312 /* of order 0 to IORDRE by U of F(Ue,Vc) or Vc=TCONST */
313 /* (fixed iso value) and Ue is the fixed extremity. */
315 /* SOMTAB: Table of NBROOT/2 sums of 2 index points */
316 /* NBROOT-II+1 and II, for II = 1, NBROOT/2. */
317 /* DIFTAB: Table of NBROOT/2 differences of 2 index points */
318 /* NBROOT-II+1 and II, for II = 1, NBROOT/2. */
320 /* OUTPUT ARGUMENTS : */
321 /* ------------------- */
322 /* SOMTAB: Table of NBROOT/2 sums of 2 index points */
323 /* NBROOT-II+1 and II, for II = 1, NBROOT/2 */
324 /* DIFTAB: Table of NBROOT/2 differences of 2 index points */
325 /* NBROOT-II+1 and II, for II = 1, NBROOT/2 */
326 /* FPNTAB: Auxiliary table. */
327 /* HERMIT: Table of coeff. 2*(IORDRE+1) Hermite polynoms */
328 /* of degree 2*IORDRE+1. */
329 /* IERCOD: Error code, */
330 /* = 0, Everythig is OK */
331 /* = 1, The value of IORDRE is out of (0,2) */
333 /* ---------------- */
335 /* REFERENCES CALLED : */
336 /* ----------------------- */
338 /* DESCRIPTION/NOTES/LIMITATIONS : */
339 /* ----------------------------------- */
340 /* The results of discretization are arranged in 2 tables */
341 /* SOMTAB and DIFTAB to earn time during the */
342 /* calculation of coefficients of the approximation curve. */
344 /* If NBROOT is uneven in SOMTAB(0,*) and DIFTAB(0,*) one stores */
345 /* the values of the median root of Legendre (0.D0 in (-1,1)). */
347 /* **********************************************************************
350 /* Name of the routine */
353 /* Parameter adjustments */
354 diftab_dim1 = *nbroot / 2 + 1;
355 diftab_offset = diftab_dim1;
356 diftab -= diftab_offset;
357 somtab_dim1 = *nbroot / 2 + 1;
358 somtab_offset = somtab_dim1;
359 somtab -= somtab_offset;
361 hermit_dim1 = (*iordre << 1) + 2;
362 hermit_offset = hermit_dim1;
363 hermit -= hermit_offset;
364 fpntab_dim1 = *nbroot;
365 fpntab_offset = fpntab_dim1 + 1;
366 fpntab -= fpntab_offset;
367 contr2_dim1 = *ndimen;
368 contr2_offset = contr2_dim1 + 1;
369 contr2 -= contr2_offset;
370 contr1_dim1 = *ndimen;
371 contr1_offset = contr1_dim1 + 1;
372 contr1 -= contr1_offset;
375 ibb = AdvApp2Var_SysBase::mnfndeb_();
377 AdvApp2Var_SysBase::mgenmsg_("MMA1CDI", 7L);
381 /* --- Recuperate 2*(IORDRE+1) coeff of 2*(IORDRE+1) of Hermite polynom ---
384 AdvApp2Var_ApproxF2var::mma1her_(iordre, &hermit[hermit_offset], iercod);
389 /* ------------------- Discretization of Hermite polynoms -----------
392 ncfhe = (*iordre + 1) << 1;
394 for (ii = 1; ii <= i__1; ++ii) {
396 for (kk = 1; kk <= i__2; ++kk) {
397 AdvApp2Var_MathBase::mmmpocur_(&ncfhe, &c__1, &ncfhe, &hermit[ii * hermit_dim1], &
398 rootlg[kk], &fpntab[kk + ii * fpntab_dim1]);
404 /* ---- Discretizations of boundary polynoms are taken ----
409 for (nd = 1; nd <= i__1; ++nd) {
411 for (ii = 1; ii <= i__2; ++ii) {
412 bid1 = contr1[nd + ii * contr1_dim1];
413 bid2 = contr2[nd + ii * contr2_dim1];
415 for (kk = 1; kk <= i__3; ++kk) {
416 kkm = nroo2 - kk + 1;
417 bid3 = bid1 * fpntab[kkm + ((ii << 1) - 1) * fpntab_dim1] +
418 bid2 * fpntab[kkm + (ii << 1) * fpntab_dim1];
419 somtab[kk + nd * somtab_dim1] -= bid3;
420 diftab[kk + nd * diftab_dim1] += bid3;
424 for (kk = 1; kk <= i__3; ++kk) {
425 kkp = (*nbroot + 1) / 2 + kk;
426 bid3 = bid1 * fpntab[kkp + ((ii << 1) - 1) * fpntab_dim1] +
427 bid2 * fpntab[kkp + (ii << 1) * fpntab_dim1];
428 somtab[kk + nd * somtab_dim1] -= bid3;
429 diftab[kk + nd * diftab_dim1] -= bid3;
437 /* ------------ Cas when discretization is done on the roots of a -----------
439 /* ---------- Legendre polynom of uneven degree, 0 is root --------
442 if (*nbroot % 2 == 1) {
444 for (nd = 1; nd <= i__1; ++nd) {
446 for (ii = 1; ii <= i__2; ++ii) {
447 bid3 = fpntab[nroo2 + 1 + ((ii << 1) - 1) * fpntab_dim1] *
448 contr1[nd + ii * contr1_dim1] + fpntab[nroo2 + 1 + (
449 ii << 1) * fpntab_dim1] * contr2[nd + ii *
453 somtab[nd * somtab_dim1] -= bid3;
454 diftab[nd * diftab_dim1] -= bid3;
461 /* ------------------------------ The End -------------------------------
463 /* --> IORDRE is not in the authorized zone. */
470 AdvApp2Var_SysBase::mgsomsg_("MMA1CDI", 7L);
475 //=======================================================================
476 //function : mma1cnt_
478 //=======================================================================
479 int mma1cnt_(integer *ndimen,
487 /* System generated locals */
488 integer contr1_dim1, contr1_offset, contr2_dim1, contr2_offset,
489 hermit_dim1, hermit_offset, crvjac_dim1, crvjac_offset, i__1,
492 /* Local variables */
493 static integer nd, ii, jj, ibb;
494 static doublereal bid;
497 /* ***********************************************************************
502 /* Add constraint to polynom. */
506 /* ALL,AB_SPECIFI::COURE&,APPROXIMATION,ADDITION,&CONSTRAINT */
508 /* INPUT ARGUMENTS : */
509 /* -------------------- */
510 /* NDIMEN: Dimension of the space */
511 /* IORDRE: Order of constraint. */
512 /* CONTR1: pt of constraint in -1, from order 0 to IORDRE. */
513 /* CONTR2: Pt of constraint in +1, from order 0 to IORDRE. */
514 /* HERMIT: Table of Hermit polynoms of order IORDRE. */
515 /* CRVJAV: Curve of approximation in Jacobi base. */
517 /* OUTPUT ARGUMENTS : */
518 /* --------------------- */
519 /* CRVJAV: Curve of approximation in Jacobi base */
520 /* to which the polynom of interpolation of constraints is added. */
523 /* ------------------ */
526 /* REFERENCES CALLED : */
527 /* --------------------- */
530 /* DESCRIPTION/NOTES/LIMITATIONS : */
531 /* ----------------------------------- */
534 /* ***********************************************************************
537 /* ***********************************************************************
539 /* Name of the routine */
541 /* ***********************************************************************
543 /* INITIALISATIONS */
544 /* ***********************************************************************
547 /* Parameter adjustments */
548 hermit_dim1 = (*iordre << 1) + 2;
549 hermit_offset = hermit_dim1;
550 hermit -= hermit_offset;
551 contr2_dim1 = *ndimen;
552 contr2_offset = contr2_dim1 + 1;
553 contr2 -= contr2_offset;
554 contr1_dim1 = *ndimen;
555 contr1_offset = contr1_dim1 + 1;
556 contr1 -= contr1_offset;
557 crvjac_dim1 = *ndgjac + 1;
558 crvjac_offset = crvjac_dim1;
559 crvjac -= crvjac_offset;
562 ibb = AdvApp2Var_SysBase::mnfndeb_();
564 AdvApp2Var_SysBase::mgenmsg_("MMA1CNT", 7L);
567 /* ***********************************************************************
570 /* ***********************************************************************
574 for (nd = 1; nd <= i__1; ++nd) {
575 i__2 = (*iordre << 1) + 1;
576 for (ii = 0; ii <= i__2; ++ii) {
579 for (jj = 1; jj <= i__3; ++jj) {
580 bid = bid + contr1[nd + jj * contr1_dim1] *
581 hermit[ii + ((jj << 1) - 1) * hermit_dim1] +
582 contr2[nd + jj * contr2_dim1] * hermit[ii + (jj << 1) * hermit_dim1];
585 crvjac[ii + nd * crvjac_dim1] = bid;
591 /* ***********************************************************************
593 /* RETURN CALLING PROGRAM */
594 /* ***********************************************************************
598 AdvApp2Var_SysBase::mgsomsg_("MMA1CNT", 7L);
604 //=======================================================================
605 //function : mma1fdi_
607 //=======================================================================
608 int mma1fdi_(integer *ndimen,
610 const AdvApp2Var_EvaluatorFunc2Var& foncnp,
624 /* System generated locals */
625 integer fpntab_dim1, somtab_dim1, somtab_offset, diftab_dim1,
626 diftab_offset, contr1_dim1, contr1_offset, contr2_dim1,
627 contr2_offset, i__1, i__2;
630 /* Local variables */
631 static integer ideb, ifin, nroo2, ideru, iderv;
632 static doublereal renor;
633 static integer ii, nd, ibb, iim, nbp, iip;
634 static doublereal bid1, bid2;
636 /* **********************************************************************
641 /* DiscretiZation of a non-polynomial function F(U,V) or of */
642 /* its derivative with fixed isoparameter. */
646 /* ALL, AB_SPECIFI::FONCTION&, DISCRETISATION, &POINT */
648 /* INPUT ARGUMENTS : */
649 /* ------------------ */
650 /* NDIMEN: Space dimension. */
651 /* UVFONC: Limits of the path of definition by U and by V of the approximated function */
652 /* FONCNP: The NAME of the non-polynomial function to be approximated */
653 /* (external program). */
654 /* ISOFAV: Fixed isoparameter for the discretization; */
655 /* = 1, discretization with fixed U and variable V. */
656 /* = 2, discretization with fixed V and variable U. */
657 /* TCONST: Iso value is also fixed. */
658 /* NBROOT: Number of INTERNAL discretization parameters. */
659 /* (if there are constraints, 2 extremities should be added).
661 /* This is also the root number of the Legendre polynom where */
662 /* the discretization is done. */
663 /* TTABLE: Table of discretization parameters and of 2 extremities */
664 /* (Respectively (-1, NBROOT Legendre roots,1) */
665 /* reframed within the adequate interval. */
666 /* IORDRE: Order of constraint imposed on the extremities of the iso. */
667 /* (If Iso-U, it is necessary to calculate the derivatives by V and vice */
669 /* = 0, the extremities of the iso are calculated. */
670 /* = 1, additionally the 1st derivative in the direction of the iso is calculated */
671 /* = 2, additionally the 2nd derivative in the direction of the iso is calculated */
672 /* IDERIV: Order of derivative transversal to fixed iso (If Iso-U=Uc */
673 /* is fixed, the derivative of order IDERIV is discretized by U of */
674 /* F(Uc,v). Same if iso-V is fixed). */
675 /* Varies from 0 (positioning) to 2 (2nd derivative). */
677 /* OUTPUT ARGUMENTS : */
678 /* ------------------- */
679 /* FPNTAB: Auxiliary table.
680 SOMTAB: Table of NBROOT/2 sums of 2 index points */
681 /* NBROOT-II+1 and II, for II = 1, NBROOT/2 */
682 /* DIFTAB: Table of NBROOT/2 differences of 2 index points */
683 /* NBROOT-II+1 and II, for II = 1, NBROOT/2 */
684 /* CONTR1: Contains, if IORDRE>=0, values IORDRE+1 in TTABLE(0)
686 /* (1st extremity) of derivatives of F(Uc,Ve) or F(Ue,Vc), */
688 /* CONTR2: Contains, if IORDRE>=0, values IORDRE+1 in */
689 /* TTABLE(NBROOT+1) (2nd extremity) of: */
690 /* If ISOFAV=1, derived of order IDERIV by U, derived */
691 /* ordre 0 to IORDRE by V of F(Uc,Ve) or Uc=TCONST */
692 /* (fixed iso value) and Ve is the fixed extremity. */
693 /* If ISOFAV=2, derivative of order IDERIV by V, derivative */
694 /* of order 0 to IORDRE by U of F(Ue,Vc) or Vc=TCONST */
695 /* (fixed iso value) and Ue is the fixed extremity. */
696 /* IERCOD: Error code > 100; Pb in evaluation of FONCNP, */
697 /* the returned error code is equal to error code of FONCNP + 100. */
700 /* ---------------- */
702 /* REFERENCES CALLED : */
703 /* ----------------------- */
705 /* DESCRIPTION/NOTES/LIMITATIONS : */
706 /* ----------------------------------- */
707 /* The results of discretization are arranged in 2 tables */
708 /* SOMTAB and DIFTAB to earn time during the */
709 /* calculation of coefficients of the approximation curve. */
711 /* If NBROOT is uneven in SOMTAB(0,*) and DIFTAB(0,*) one stores */
712 /* the values of the median root of Legendre (0.D0 in (-1,1)). */
714 /* Function F(u,v) defined in UVFONC is reparameterized in */
715 /* (-1,1)x(-1,1). Then 1st and 2nd derivatives are renormalized. */
718 /* **********************************************************************
721 /* Name of the routine */
724 /* Parameter adjustments */
726 diftab_dim1 = *nbroot / 2 + 1;
727 diftab_offset = diftab_dim1;
728 diftab -= diftab_offset;
729 somtab_dim1 = *nbroot / 2 + 1;
730 somtab_offset = somtab_dim1;
731 somtab -= somtab_offset;
732 fpntab_dim1 = *ndimen;
734 contr2_dim1 = *ndimen;
735 contr2_offset = contr2_dim1 + 1;
736 contr2 -= contr2_offset;
737 contr1_dim1 = *ndimen;
738 contr1_offset = contr1_dim1 + 1;
739 contr1 -= contr1_offset;
742 ibb = AdvApp2Var_SysBase::mnfndeb_();
744 AdvApp2Var_SysBase::mgenmsg_("MMA1FDI", 7L);
748 /* --------------- Definition of the nb of points to calculate --------------
750 /* --> If constraints, the limits are also taken */
754 /* --> Otherwise, only Legendre roots (reframed) are used
760 /* --> Nb of point to calculate. */
761 nbp = ifin - ideb + 1;
764 /* --------------- Determination of the order of global derivation --------
766 /* --> ISOFAV takes only values 1 or 2. */
767 /* if Iso-U, derive by U of order IDERIV */
771 d__1 = (uvfonc[4] - uvfonc[3]) / 2.;
772 renor = AdvApp2Var_MathBase::pow__di(&d__1, ideriv);
773 /* if Iso-V, derive by V of order IDERIV */
777 d__1 = (uvfonc[6] - uvfonc[5]) / 2.;
778 renor = AdvApp2Var_MathBase::pow__di(&d__1, ideriv);
781 /* ----------- Discretization on roots of the ---------------
783 /* ---------------------- Legendre polynom of degree NBROOT -------------------
786 foncnp.Evaluate (ndimen,
795 &fpntab[ideb * fpntab_dim1 + 1],
801 for (nd = 1; nd <= i__1; ++nd) {
803 for (ii = 1; ii <= i__2; ++ii) {
804 iip = (*nbroot + 1) / 2 + ii;
805 iim = nroo2 - ii + 1;
806 bid1 = fpntab[nd + iim * fpntab_dim1];
807 bid2 = fpntab[nd + iip * fpntab_dim1];
808 somtab[ii + nd * somtab_dim1] = renor * (bid2 + bid1);
809 diftab[ii + nd * diftab_dim1] = renor * (bid2 - bid1);
815 /* ------------ Case when discretisation is done on roots of a ----
817 /* ---------- Legendre polynom of uneven degree, 0 is root --------
820 if (*nbroot % 2 == 1) {
822 for (nd = 1; nd <= i__1; ++nd) {
823 somtab[nd * somtab_dim1] = renor * fpntab[nd + (nroo2 + 1) *
825 diftab[nd * diftab_dim1] = renor * fpntab[nd + (nroo2 + 1) *
831 for (nd = 1; nd <= i__1; ++nd) {
832 somtab[nd * somtab_dim1] = 0.;
833 diftab[nd * diftab_dim1] = 0.;
838 /* --------------------- Take into account constraints ----------------
842 /* --> Recover already calculated extremities. */
844 for (nd = 1; nd <= i__1; ++nd) {
845 contr1[nd + contr1_dim1] = renor * fpntab[nd];
846 contr2[nd + contr2_dim1] = renor * fpntab[nd + (*nbroot + 1) *
850 /* --> Nb of points to calculate/call to FONCNP */
852 /* If Iso-U, derive by V till order IORDRE */
854 /* --> Factor of normalisation 1st derivative. */
855 bid1 = (uvfonc[6] - uvfonc[5]) / 2.;
857 for (iderv = 1; iderv <= i__1; ++iderv) {
858 foncnp.Evaluate (ndimen, &uvfonc[3], &uvfonc[5], isofav, tconst,
859 &nbp, ttable, &ideru, &iderv, &contr1[(iderv + 1) *
860 contr1_dim1 + 1], iercod);
867 for (iderv = 1; iderv <= i__1; ++iderv) {
868 foncnp.Evaluate (ndimen, &uvfonc[3], &uvfonc[5], isofav, tconst,
869 &nbp, &ttable[*nbroot + 1], &ideru, &iderv, &contr2[(
870 iderv + 1) * contr2_dim1 + 1], iercod);
876 /* If Iso-V, derive by U till order IORDRE */
878 /* --> Factor of normalization 1st derivative. */
879 bid1 = (uvfonc[4] - uvfonc[3]) / 2.;
881 for (ideru = 1; ideru <= i__1; ++ideru) {
882 foncnp.Evaluate (ndimen, &uvfonc[3], &uvfonc[5], isofav, tconst,
883 &nbp, ttable, &ideru, &iderv, &contr1[(ideru + 1) *
884 contr1_dim1 + 1], iercod);
891 for (ideru = 1; ideru <= i__1; ++ideru) {
892 foncnp.Evaluate (ndimen, &uvfonc[3], &uvfonc[5], isofav, tconst,
893 &nbp, &ttable[*nbroot + 1], &ideru, &iderv, &contr2[(
894 ideru + 1) * contr2_dim1 + 1], iercod);
902 /* ------------------------- Normalization of derivatives -------------
904 /* (The function is redefined on (-1,1)*(-1,1)) */
907 for (ii = 1; ii <= i__1; ++ii) {
910 for (nd = 1; nd <= i__2; ++nd) {
911 contr1[nd + (ii + 1) * contr1_dim1] *= bid2;
912 contr2[nd + (ii + 1) * contr2_dim1] *= bid2;
919 /* ------------------------------ The end -------------------------------
925 AdvApp2Var_SysBase::maermsg_("MMA1FDI", iercod, 7L);
928 AdvApp2Var_SysBase::mgsomsg_("MMA1FDI", 7L);
933 //=======================================================================
934 //function : mma1fer_
936 //=======================================================================
937 int mma1fer_(integer *,//ndimen,
951 /* System generated locals */
952 integer crvjac_dim1, crvjac_offset, i__1, i__2;
954 /* Local variables */
955 static integer idim, ncfja, ncfnw, ndses, ii, kk, ibb, ier;
959 /* ***********************************************************************
964 /* Calculate the degree and the errors of approximation of a border. */
968 /* TOUS,AB_SPECIFI :: COURBE&,TRONCATURE, &PRECISION */
970 /* INPUT ARGUMENTS : */
971 /* -------------------- */
973 /* NDIMEN: Total Dimension of the space (sum of dimensions of sub-spaces) */
974 /* NBSESP: Number of "independent" sub-spaces. */
975 /* NDIMSE: Table of dimensions of sub-spaces. */
976 /* IORDRE: Order of constraint at the extremities of the border */
977 /* -1 = no constraints, */
978 /* 0 = constraints of passage to limits (i.e. C0), */
979 /* 1 = C0 + constraintes of 1st derivatives (i.e. C1), */
980 /* 2 = C1 + constraintes of 2nd derivatives (i.e. C2). */
981 /* NDGJAC: Degree of development in series to use for the calculation
982 /* in the base of Jacobi. */
983 /* CRVJAC: Table of coeff. of the curve of approximation in the */
984 /* base of Jacobi. */
985 /* NCFLIM: Max number of coeff of the polynomial curve */
986 /* of approximation (should be above or equal to */
987 /* 2*IORDRE+2 and below or equal to 50). */
988 /* EPSAPR: Table of errors of approximations that cannot be passed, */
989 /* sub-space by sub-space. */
991 /* OUTPUT ARGUMENTS : */
992 /* --------------------- */
993 /* YCVMAX: Auxiliary Table. */
994 /* ERRMAX: Table of errors (sub-space by sub-space) */
995 /* MAXIMUM made in the approximation of FONCNP by */
997 /* ERRMOY: Table of errors (sub-space by sub-space) */
998 /* AVERAGE made in the approximation of FONCNP by */
1000 /* NCOEFF: Number of significative coeffs. of the calculated "curve". */
1001 /* IERCOD: Error code */
1003 /* =-1, warning, required tolerance can't be */
1004 /* met with coefficients NFCLIM. */
1005 /* = 1, order of constraints (IORDRE) is not within authorised values */
1008 /* COMMONS USED : */
1009 /* ------------------ */
1011 /* REFERENCES CALLED : */
1012 /* --------------------- */
1014 /* DESCRIPTION/NOTES/LIMITATIONS : */
1015 /* ----------------------------------- */
1017 /* **********************************************************************
1020 /* Name of the routine */
1023 /* Parameter adjustments */
1029 crvjac_dim1 = *ndgjac + 1;
1030 crvjac_offset = crvjac_dim1;
1031 crvjac -= crvjac_offset;
1034 ibb = AdvApp2Var_SysBase::mnfndeb_();
1036 AdvApp2Var_SysBase::mgenmsg_("MMA1FER", 7L);
1041 ncfja = *ndgjac + 1;
1043 /* ------------ Calculate the degree of the curve and of the Max error --------
1045 /* -------------- of approximation for all sub-spaces --------
1049 for (ii = 1; ii <= i__1; ++ii) {
1052 /* ------------ cutting of coeff. and calculation of Max error -------
1055 AdvApp2Var_MathBase::mmtrpjj_(&ncfja, &ndses, &ncfja, &epsapr[ii], iordre, &crvjac[idim *
1056 crvjac_dim1], &ycvmax[1], &errmax[ii], &ncfnw);
1058 /* ******************************************************************
1060 /* ------------- If precision OK, calculate the average error -------
1062 /* ******************************************************************
1065 if (ncfnw <= *ncflim) {
1066 mmaperm_(&ncfja, &ndses, &ncfja, iordre, &crvjac[idim *
1067 crvjac_dim1], &ncfnw, &errmoy[ii]);
1068 *ncoeff = advapp_max(ncfnw,*ncoeff);
1070 /* ------------- Set the declined coefficients to 0.D0 -----------
1073 nbr0 = *ncflim - ncfnw;
1076 for (kk = 1; kk <= i__2; ++kk) {
1077 AdvApp2Var_SysBase::mvriraz_(&nbr0,
1078 (char *)&crvjac[ncfnw + (idim + kk - 1) * crvjac_dim1]);
1084 /* **************************************************************
1086 /* ------------------- If required precision can't be reached----
1088 /* **************************************************************
1093 /* ------------------------- calculate the Max error ------------
1096 AdvApp2Var_MathBase::mmaperx_(&ncfja, &ndses, &ncfja, iordre, &crvjac[idim *
1097 crvjac_dim1], ncflim, &ycvmax[1], &errmax[ii], &ier);
1102 /* -------------------- nb of coeff to be returned -------------
1107 /* ------------------- and calculation of the average error ----
1110 mmaperm_(&ncfja, &ndses, &ncfja, iordre, &crvjac[idim *
1111 crvjac_dim1], ncflim, &errmoy[ii]);
1119 /* ------------------------------ The end -------------------------------
1121 /* --> The order of constraints is not within autorized values. */
1128 AdvApp2Var_SysBase::maermsg_("MMA1FER", iercod, 7L);
1131 AdvApp2Var_SysBase::mgsomsg_("MMA1FER", 7L);
1137 //=======================================================================
1138 //function : mma1her_
1140 //=======================================================================
1141 int AdvApp2Var_ApproxF2var::mma1her_(const integer *iordre,
1145 /* System generated locals */
1146 integer hermit_dim1, hermit_offset;
1148 /* Local variables */
1153 /* **********************************************************************
1158 /* Calculate 2*(IORDRE+1) Hermit polynoms of degree 2*IORDRE+1 */
1163 /* ALL, AB_SPECIFI::CONTRAINTE&, INTERPOLATION, &POLYNOME */
1165 /* INPUT ARGUMENTS : */
1166 /* ------------------ */
1167 /* IORDRE: Order of constraint. */
1168 /* = 0, Polynom of interpolation of order C0 on (-1,1). */
1169 /* = 1, Polynom of interpolation of order C0 and C1 on (-1,1). */
1170 /* = 2, Polynom of interpolation of order C0, C1 and C2 on (-1,1).
1173 /* OUTPUT ARGUMENTS : */
1174 /* ------------------- */
1175 /* HERMIT: Table of 2*IORDRE+2 coeff. of each of 2*(IORDRE+1) */
1176 /* HERMIT polynom. */
1177 /* IERCOD: Error code, */
1179 /* = 1, required order of constraint is not managed here. */
1180 /* COMMONS USED : */
1181 /* ---------------- */
1183 /* REFERENCES CALLED : */
1184 /* ----------------------- */
1186 /* DESCRIPTION/NOTES/LIMITATIONS : */
1187 /* ----------------------------------- */
1188 /* The part of HERMIT(*,2*i+j) table where j=1 or 2 and i=0 to IORDRE,
1189 /* contains the coefficients of the polynom of degree 2*IORDRE+1 */
1190 /* such as ALL values in -1 and in +1 of this polynom and its */
1191 /* derivatives till order of derivation IORDRE are NULL, */
1192 /* EXCEPT for the derivative of order i: */
1193 /* - valued 1 in -1 if j=1 */
1194 /* - valued 1 in +1 if j=2. */
1196 /* **********************************************************************
1199 /* Name of the routine */
1202 /* Parameter adjustments */
1203 hermit_dim1 = (*iordre + 1) << 1;
1204 hermit_offset = hermit_dim1 + 1;
1205 hermit -= hermit_offset;
1208 ibb = AdvApp2Var_SysBase::mnfndeb_();
1210 AdvApp2Var_SysBase::mgenmsg_("MMA1HER", 7L);
1214 /* --- Recover (IORDRE+2) coeff of 2*(IORDRE+1) Hermit polynoms --
1218 hermit[hermit_dim1 + 1] = .5;
1219 hermit[hermit_dim1 + 2] = -.5;
1221 hermit[(hermit_dim1 << 1) + 1] = .5;
1222 hermit[(hermit_dim1 << 1) + 2] = .5;
1223 } else if (*iordre == 1) {
1224 hermit[hermit_dim1 + 1] = .5;
1225 hermit[hermit_dim1 + 2] = -.75;
1226 hermit[hermit_dim1 + 3] = 0.;
1227 hermit[hermit_dim1 + 4] = .25;
1229 hermit[(hermit_dim1 << 1) + 1] = .5;
1230 hermit[(hermit_dim1 << 1) + 2] = .75;
1231 hermit[(hermit_dim1 << 1) + 3] = 0.;
1232 hermit[(hermit_dim1 << 1) + 4] = -.25;
1234 hermit[hermit_dim1 * 3 + 1] = .25;
1235 hermit[hermit_dim1 * 3 + 2] = -.25;
1236 hermit[hermit_dim1 * 3 + 3] = -.25;
1237 hermit[hermit_dim1 * 3 + 4] = .25;
1239 hermit[(hermit_dim1 << 2) + 1] = -.25;
1240 hermit[(hermit_dim1 << 2) + 2] = -.25;
1241 hermit[(hermit_dim1 << 2) + 3] = .25;
1242 hermit[(hermit_dim1 << 2) + 4] = .25;
1243 } else if (*iordre == 2) {
1244 hermit[hermit_dim1 + 1] = .5;
1245 hermit[hermit_dim1 + 2] = -.9375;
1246 hermit[hermit_dim1 + 3] = 0.;
1247 hermit[hermit_dim1 + 4] = .625;
1248 hermit[hermit_dim1 + 5] = 0.;
1249 hermit[hermit_dim1 + 6] = -.1875;
1251 hermit[(hermit_dim1 << 1) + 1] = .5;
1252 hermit[(hermit_dim1 << 1) + 2] = .9375;
1253 hermit[(hermit_dim1 << 1) + 3] = 0.;
1254 hermit[(hermit_dim1 << 1) + 4] = -.625;
1255 hermit[(hermit_dim1 << 1) + 5] = 0.;
1256 hermit[(hermit_dim1 << 1) + 6] = .1875;
1258 hermit[hermit_dim1 * 3 + 1] = .3125;
1259 hermit[hermit_dim1 * 3 + 2] = -.4375;
1260 hermit[hermit_dim1 * 3 + 3] = -.375;
1261 hermit[hermit_dim1 * 3 + 4] = .625;
1262 hermit[hermit_dim1 * 3 + 5] = .0625;
1263 hermit[hermit_dim1 * 3 + 6] = -.1875;
1265 hermit[(hermit_dim1 << 2) + 1] = -.3125;
1266 hermit[(hermit_dim1 << 2) + 2] = -.4375;
1267 hermit[(hermit_dim1 << 2) + 3] = .375;
1268 hermit[(hermit_dim1 << 2) + 4] = .625;
1269 hermit[(hermit_dim1 << 2) + 5] = -.0625;
1270 hermit[(hermit_dim1 << 2) + 6] = -.1875;
1272 hermit[hermit_dim1 * 5 + 1] = .0625;
1273 hermit[hermit_dim1 * 5 + 2] = -.0625;
1274 hermit[hermit_dim1 * 5 + 3] = -.125;
1275 hermit[hermit_dim1 * 5 + 4] = .125;
1276 hermit[hermit_dim1 * 5 + 5] = .0625;
1277 hermit[hermit_dim1 * 5 + 6] = -.0625;
1279 hermit[hermit_dim1 * 6 + 1] = .0625;
1280 hermit[hermit_dim1 * 6 + 2] = .0625;
1281 hermit[hermit_dim1 * 6 + 3] = -.125;
1282 hermit[hermit_dim1 * 6 + 4] = -.125;
1283 hermit[hermit_dim1 * 6 + 5] = .0625;
1284 hermit[hermit_dim1 * 6 + 6] = .0625;
1289 /* ------------------------------ The End -------------------------------
1292 AdvApp2Var_SysBase::maermsg_("MMA1HER", iercod, 7L);
1294 AdvApp2Var_SysBase::mgsomsg_("MMA1HER", 7L);
1298 //=======================================================================
1299 //function : mma1jak_
1301 //=======================================================================
1302 int mma1jak_(integer *ndimen,
1312 /* System generated locals */
1313 integer somtab_dim1, somtab_offset, diftab_dim1, diftab_offset,
1314 crvjac_dim1, crvjac_offset, cgauss_dim1;
1316 /* Local variables */
1319 /* **********************************************************************
1324 /* Calculate the curve of approximation of a non-polynomial function */
1325 /* in the base of Jacobi. */
1329 /* FUNCTION,DISCRETISATION,APPROXIMATION,CONSTRAINT,CURVE,JACOBI */
1331 /* INPUT ARGUMENTS : */
1332 /* ------------------ */
1333 /* NDIMEN: Total dimension of the space (sum of dimensions */
1334 /* of sub-spaces) */
1335 /* NBROOT: Nb of points of discretization of the iso, extremities not
1337 /* IORDRE: Order of constraint at the extremities of the boundary */
1338 /* -1 = no constraints, */
1339 /* 0 = constraints of passage of limits (i.e. C0), */
1340 /* 1 = C0 + constraints of 1st derivatives (i.e. C1), */
1341 /* 2 = C1 + constraints of 2nd derivatives (i.e. C2). */
1342 /* NDGJAC: Degree of development in series to be used for calculation in the
1343 /* base of Jacobi. */
1345 /* OUTPUT ARGUMENTS : */
1346 /* ------------------- */
1347 /* CRVJAC : Curve of approximation of FONCNP with (eventually) */
1348 /* taking into account of constraints at the extremities. */
1349 /* This curve is of degree NDGJAC. */
1350 /* IERCOD : Error code : */
1351 /* 0 = All is ok. */
1352 /* 33 = Pb to return data of du block data */
1353 /* of coeff. of integration by GAUSS method. */
1354 /* by program MMAPPTT. */
1356 /* COMMONS USED : */
1357 /* ---------------- */
1359 /* REFERENCES CALLED : */
1360 /* ----------------------- */
1362 /* **********************************************************************
1365 /* Name of the routine */
1367 /* Parameter adjustments */
1368 diftab_dim1 = *nbroot / 2 + 1;
1369 diftab_offset = diftab_dim1;
1370 diftab -= diftab_offset;
1371 somtab_dim1 = *nbroot / 2 + 1;
1372 somtab_offset = somtab_dim1;
1373 somtab -= somtab_offset;
1374 crvjac_dim1 = *ndgjac + 1;
1375 crvjac_offset = crvjac_dim1;
1376 crvjac -= crvjac_offset;
1377 cgauss_dim1 = *nbroot / 2 + 1;
1380 ibb = AdvApp2Var_SysBase::mnfndeb_();
1382 AdvApp2Var_SysBase::mgenmsg_("MMA1JAK", 7L);
1386 /* ----------------- Recover coeffs of integration by Gauss -----------
1389 AdvApp2Var_ApproxF2var::mmapptt_(ndgjac, nbroot, iordre, cgauss, iercod);
1395 /* --------------- Calculate the curve in the base of Jacobi -----------
1398 mmmapcoe_(ndimen, ndgjac, iordre, nbroot, &somtab[somtab_offset], &diftab[
1399 diftab_offset], cgauss, &crvjac[crvjac_offset]);
1401 /* ------------------------------ The End -------------------------------
1406 AdvApp2Var_SysBase::maermsg_("MMA1JAK", iercod, 7L);
1409 AdvApp2Var_SysBase::mgsomsg_("MMA1JAK", 7L);
1414 //=======================================================================
1415 //function : mma1noc_
1417 //=======================================================================
1418 int mma1noc_(doublereal *dfuvin,
1427 /* System generated locals */
1431 /* Local variables */
1432 static doublereal rider, riord;
1433 static integer nd, ibb;
1434 static doublereal bid;
1435 /* **********************************************************************
1440 /* Normalization of constraints of derivatives, defined on DFUVIN */
1441 /* on block DUVOUT. */
1445 /* ALL, AB_SPECIFI::VECTEUR&,DERIVEE&,NORMALISATION,&VECTEUR */
1447 /* INPUT ARGUMENTS : */
1448 /* ------------------ */
1449 /* DFUVIN: Limits of the block of definition by U and by V where
1451 /* constraints CNTRIN are defined. */
1452 /* NDIMEN: Dimension of the space. */
1453 /* IORDRE: Order of constraint imposed at the extremities of the iso. */
1454 /* (if Iso-U, it is necessary to calculate derivatives by V and vice */
1456 /* = 0, the extremities of the iso are calculated */
1457 /* = 1, additionally the 1st derivative in the direction */
1458 /* of the iso is calculated */
1459 /* = 2, additionally the 2nd derivative in the direction */
1460 /* of the iso is calculated */
1461 /* CNTRIN: Contains, if IORDRE>=0, IORDRE+1 derivatives */
1462 /* of order IORDRE of F(Uc,v) or of F(u,Vc), following the */
1463 /* value of ISOFAV, RENORMALIZED by u and v in (-1,1). */
1464 /* DUVOUT: Limits of the block of definition by U and by V where the */
1465 /* constraints CNTOUT will be defined. */
1466 /* ISOFAV: Isoparameter fixed for the discretization; */
1467 /* = 1, discretization with fixed U=Uc and variable V. */
1468 /* = 2, discretization with fixed V=Vc and variable U. */
1469 /* IDERIV: Ordre de derivee transverse a l'iso fixee (Si Iso-U=Uc */
1470 /* is fixed, the derivative of order IDERIV is discretized by U */
1471 /* of F(Uc,v). The same if iso-V is fixed). */
1472 /* Varies from (positioning) to 2 (2nd derivative). */
1474 /* OUTPUT ARGUMENTS : */
1475 /* ------------------- */
1476 /* CNTOUT: Contains, if IORDRE>=0, IORDRE+1 derivatives */
1477 /* of order IORDRE of F(Uc,v) or of F(u,Vc), depending on the */
1478 /* value of ISOFAV, RENORMALIZED for u and v in DUVOUT. */
1480 /* COMMONS USED : */
1481 /* ---------------- */
1483 /* REFERENCES CALLED : */
1484 /* --------------------- */
1486 /* DESCRIPTION/NOTES/LIMITATIONS : */
1487 /* ------------------------------- */
1488 /* CNTRIN can be an output/input argument, */
1491 /* CALL MMA1NOC(DFUVIN,NDIMEN,IORDRE,CNTRIN,DUVOUT */
1492 /* 1 ,ISOFAV,IDERIV,CNTRIN) */
1496 /* **********************************************************************
1499 /* Name of the routine */
1502 /* Parameter adjustments */
1509 ibb = AdvApp2Var_SysBase::mnfndeb_();
1511 AdvApp2Var_SysBase::mgenmsg_("MMA1NOC", 7L);
1514 /* --------------- Determination of coefficients of normalization -------
1518 d__1 = (dfuvin[4] - dfuvin[3]) / (duvout[4] - duvout[3]);
1519 rider = AdvApp2Var_MathBase::pow__di(&d__1, ideriv);
1520 d__1 = (dfuvin[6] - dfuvin[5]) / (duvout[6] - duvout[5]);
1521 riord = AdvApp2Var_MathBase::pow__di(&d__1, iordre);
1524 d__1 = (dfuvin[6] - dfuvin[5]) / (duvout[6] - duvout[5]);
1525 rider = AdvApp2Var_MathBase::pow__di(&d__1, ideriv);
1526 d__1 = (dfuvin[4] - dfuvin[3]) / (duvout[4] - duvout[3]);
1527 riord = AdvApp2Var_MathBase::pow__di(&d__1, iordre);
1530 /* ------------- Renormalization of the vector of constraint ---------------
1533 bid = rider * riord;
1535 for (nd = 1; nd <= i__1; ++nd) {
1536 cntout[nd] = bid * cntrin[nd];
1540 /* ------------------------------ The end -------------------------------
1544 AdvApp2Var_SysBase::mgsomsg_("MMA1NOC", 7L);
1549 //=======================================================================
1550 //function : mma1nop_
1552 //=======================================================================
1553 int mma1nop_(integer *nbroot,
1561 /* System generated locals */
1564 /* Local variables */
1565 static doublereal alinu, blinu, alinv, blinv;
1566 static integer ii, ibb;
1568 /* ***********************************************************************
1573 /* Normalization of parameters of an iso, starting from */
1574 /* parametric block and parameters on (-1,1). */
1578 /* TOUS,AB_SPECIFI :: ISO&,POINT&,NORMALISATION,&POINT,&ISO */
1580 /* INPUT ARGUMENTS : */
1581 /* -------------------- */
1582 /* NBROOT: Nb of points of discretisation INSIDE the iso */
1583 /* defined on (-1,1). */
1584 /* ROOTLG: Table of discretization parameters on )-1,1( */
1586 /* UVFONC: Block of definition of the iso */
1587 /* ISOFAV: = 1, this is iso-u; =2, this is iso-v. */
1589 /* OUTPUT ARGUMENTS : */
1590 /* --------------------- */
1591 /* TTABLE: Table of parameters renormalized on UVFONC of the iso.
1593 /* IERCOD: = 0, OK */
1594 /* = 1, ISOFAV is out of allowed values. */
1597 /* **********************************************************************
1599 /* Name of the routine */
1602 /* Parameter adjustments */
1607 ibb = AdvApp2Var_SysBase::mnfndeb_();
1609 AdvApp2Var_SysBase::mgenmsg_("MMA1NOP", 7L);
1612 alinu = (uvfonc[4] - uvfonc[3]) / 2.;
1613 blinu = (uvfonc[4] + uvfonc[3]) / 2.;
1614 alinv = (uvfonc[6] - uvfonc[5]) / 2.;
1615 blinv = (uvfonc[6] + uvfonc[5]) / 2.;
1618 ttable[0] = uvfonc[5];
1620 for (ii = 1; ii <= i__1; ++ii) {
1621 ttable[ii] = alinv * rootlg[ii] + blinv;
1624 ttable[*nbroot + 1] = uvfonc[6];
1625 } else if (*isofav == 2) {
1626 ttable[0] = uvfonc[3];
1628 for (ii = 1; ii <= i__1; ++ii) {
1629 ttable[ii] = alinu * rootlg[ii] + blinu;
1632 ttable[*nbroot + 1] = uvfonc[4];
1639 /* ------------------------------ THE END -------------------------------
1648 AdvApp2Var_SysBase::maermsg_("MMA1NOP", iercod, 7L);
1651 AdvApp2Var_SysBase::mgsomsg_("MMA1NOP", 7L);
1658 //=======================================================================
1659 //function : mma2ac1_
1661 //=======================================================================
1662 int AdvApp2Var_ApproxF2var::mma2ac1_(integer const *ndimen,
1663 integer const *mxujac,
1664 integer const *mxvjac,
1665 integer const *iordru,
1666 integer const *iordrv,
1667 doublereal const *contr1,
1668 doublereal const * contr2,
1669 doublereal const *contr3,
1670 doublereal const *contr4,
1671 doublereal const *uhermt,
1672 doublereal const *vhermt,
1676 /* System generated locals */
1677 integer contr1_dim1, contr1_dim2, contr1_offset, contr2_dim1, contr2_dim2,
1678 contr2_offset, contr3_dim1, contr3_dim2, contr3_offset,
1679 contr4_dim1, contr4_dim2, contr4_offset, uhermt_dim1,
1680 uhermt_offset, vhermt_dim1, vhermt_offset, patjac_dim1,
1681 patjac_dim2, patjac_offset, i__1, i__2, i__3, i__4, i__5;
1683 /* Local variables */
1684 static logical ldbg;
1685 static integer ndgu, ndgv;
1686 static doublereal bidu1, bidu2, bidv1, bidv2;
1687 static integer ioru1, iorv1, ii, nd, jj, ku, kv;
1688 static doublereal cnt1, cnt2, cnt3, cnt4;
1690 /* **********************************************************************
1695 /* Add polynoms of edge constraints. */
1699 /* TOUS,AB_SPECIFI::POINT&,CONTRAINTE&,ADDITION,&POLYNOME */
1701 /* INPUT ARGUMENTS : */
1702 /* ------------------ */
1703 /* NDIMEN: Dimension of the space. */
1704 /* MXUJAC: Max degree of the polynom of approximation by U. The */
1705 /* representation in the orthogonal base starts from degree */
1706 /* 0 to degree MXUJAC-2*(IORDRU+1). The polynomial base is the */
1707 /* base of Jacobi of order -1 (Legendre), 0, 1 or 2 */
1708 /* MXVJAC: Max degree of the polynom of approximation by V. The */
1709 /* representation in the orthogonal base starts from degree */
1710 /* 0 to degree MXUJAC-2*(IORDRU+1). The polynomial base is the */
1711 /* base of Jacobi of order -1 (Legendre), 0, 1 or 2 */
1712 /* IORDRU: Order of the base of Jacobi (-1,0,1 or 2) by U. Corresponds */
1713 /* to the step of constraints: C0, C1 or C2. */
1714 /* IORDRV: Order of the base of Jacobi (-1,0,1 or 2) by V. Corresponds */
1715 /* to the step of constraints: C0, C1 or C2. */
1716 /* CONTR1: Contains, if IORDRU and IORDRV>=0, the values at the */
1717 /* extremities of F(U0,V0) and its derivatives. */
1718 /* CONTR2: Contains, if IORDRU and IORDRV>=0, the values at the */
1719 /* extremities of F(U1,V0) and its derivatives. */
1720 /* CONTR3: Contains, if IORDRU and IORDRV>=0, the values at the */
1721 /* extremities of F(U0,V1) and its derivatives. */
1722 /* CONTR4: Contains, if IORDRU and IORDRV>=0, the values at the */
1723 /* extremities of F(U1,V1) and its derivatives. */
1724 /* UHERMT: Coeff. of Hermit polynoms of order IORDRU. */
1725 /* VHERMT: Coeff. of Hermit polynoms of order IORDRV. */
1726 /* PATJAC: Table of coefficients of the polynom P(u,v) of approximation */
1727 /* of F(u,v) WITHOUT taking into account the constraints. */
1729 /* OUTPUT ARGUMENTS : */
1730 /* ------------------- */
1731 /* PATJAC: Table of coefficients of the polynom P(u,v) by approximation */
1732 /* of F(u,v) WITH taking into account of constraints. */
1734 /* **********************************************************************
1736 /* Name of the routine */
1738 /* --------------------------- Initialization --------------------------
1741 /* Parameter adjustments */
1742 patjac_dim1 = *mxujac + 1;
1743 patjac_dim2 = *mxvjac + 1;
1744 patjac_offset = patjac_dim1 * patjac_dim2;
1745 patjac -= patjac_offset;
1746 uhermt_dim1 = (*iordru << 1) + 2;
1747 uhermt_offset = uhermt_dim1;
1748 uhermt -= uhermt_offset;
1749 vhermt_dim1 = (*iordrv << 1) + 2;
1750 vhermt_offset = vhermt_dim1;
1751 vhermt -= vhermt_offset;
1752 contr4_dim1 = *ndimen;
1753 contr4_dim2 = *iordru + 2;
1754 contr4_offset = contr4_dim1 * (contr4_dim2 + 1) + 1;
1755 contr4 -= contr4_offset;
1756 contr3_dim1 = *ndimen;
1757 contr3_dim2 = *iordru + 2;
1758 contr3_offset = contr3_dim1 * (contr3_dim2 + 1) + 1;
1759 contr3 -= contr3_offset;
1760 contr2_dim1 = *ndimen;
1761 contr2_dim2 = *iordru + 2;
1762 contr2_offset = contr2_dim1 * (contr2_dim2 + 1) + 1;
1763 contr2 -= contr2_offset;
1764 contr1_dim1 = *ndimen;
1765 contr1_dim2 = *iordru + 2;
1766 contr1_offset = contr1_dim1 * (contr1_dim2 + 1) + 1;
1767 contr1 -= contr1_offset;
1770 ldbg = AdvApp2Var_SysBase::mnfndeb_() >= 3;
1772 AdvApp2Var_SysBase::mgenmsg_("MMA2AC1", 7L);
1775 /* ------------ SUBTRACTION OF ANGULAR CONSTRAINTS -------------------
1778 ioru1 = *iordru + 1;
1779 iorv1 = *iordrv + 1;
1780 ndgu = (*iordru << 1) + 1;
1781 ndgv = (*iordrv << 1) + 1;
1784 for (jj = 1; jj <= i__1; ++jj) {
1786 for (ii = 1; ii <= i__2; ++ii) {
1788 for (nd = 1; nd <= i__3; ++nd) {
1789 cnt1 = contr1[nd + (ii + jj * contr1_dim2) * contr1_dim1];
1790 cnt2 = contr2[nd + (ii + jj * contr2_dim2) * contr2_dim1];
1791 cnt3 = contr3[nd + (ii + jj * contr3_dim2) * contr3_dim1];
1792 cnt4 = contr4[nd + (ii + jj * contr4_dim2) * contr4_dim1];
1794 for (kv = 0; kv <= i__4; ++kv) {
1795 bidv1 = vhermt[kv + ((jj << 1) - 1) * vhermt_dim1];
1796 bidv2 = vhermt[kv + (jj << 1) * vhermt_dim1];
1798 for (ku = 0; ku <= i__5; ++ku) {
1799 bidu1 = uhermt[ku + ((ii << 1) - 1) * uhermt_dim1];
1800 bidu2 = uhermt[ku + (ii << 1) * uhermt_dim1];
1801 patjac[ku + (kv + nd * patjac_dim2) * patjac_dim1] =
1802 patjac[ku + (kv + nd * patjac_dim2) *
1803 patjac_dim1] - bidu1 * bidv1 * cnt1 - bidu2 *
1804 bidv1 * cnt2 - bidu1 * bidv2 * cnt3 - bidu2 *
1817 /* ------------------------------ The end -------------------------------
1821 AdvApp2Var_SysBase::mgsomsg_("MMA2AC1", 7L);
1826 //=======================================================================
1827 //function : mma2ac2_
1829 //=======================================================================
1830 int AdvApp2Var_ApproxF2var::mma2ac2_(const integer *ndimen,
1831 const integer *mxujac,
1832 const integer *mxvjac,
1833 const integer *iordrv,
1834 const integer *nclimu,
1835 const integer *ncfiv1,
1836 const doublereal *crbiv1,
1837 const integer *ncfiv2,
1838 const doublereal *crbiv2,
1839 const doublereal *vhermt,
1843 /* System generated locals */
1844 integer crbiv1_dim1, crbiv1_dim2, crbiv1_offset, crbiv2_dim1, crbiv2_dim2,
1845 crbiv2_offset, patjac_dim1, patjac_dim2, patjac_offset,
1846 vhermt_dim1, vhermt_offset, i__1, i__2, i__3, i__4;
1848 /* Local variables */
1849 static logical ldbg;
1850 static integer ndgv1, ndgv2, ii, jj, nd, kk;
1851 static doublereal bid1, bid2;
1853 /* **********************************************************************
1858 /* Add polynoms of constraints */
1862 /* FUNCTION,APPROXIMATION,COEFFICIENT,POLYNOM */
1864 /* INPUT ARGUMENTS : */
1865 /* ------------------ */
1866 /* NDIMEN: Dimension of the space. */
1867 /* MXUJAC: Max degree of the polynom of approximation by U. The */
1868 /* representation in the orthogonal base starts from degree */
1869 /* 0 to degree MXUJAC-2*(IORDRU+1). The polynomial base is the */
1870 /* base of Jacobi of order -1 (Legendre), 0, 1 or 2 */
1871 /* MXVJAC: Max degree of the polynom of approximation by V. The */
1872 /* representation in the orthogonal base starts from degree */
1873 /* 0 to degree MXUJAC-2*(IORDRU+1). The polynomial base is the */
1874 /* base of Jacobi of order -1 (Legendre), 0, 1 or 2 */
1875 /* IORDRV: Order of the base of Jacobi (-1,0,1 or 2) by V. Corresponds */
1876 /* to the step of constraints: C0, C1 or C2. */
1877 /* NCLIMU LIMIT nb of coeff by u of the solution P(u,v)
1878 * NCFIV1: Nb of Coeff. of curves stored in CRBIV1. */
1879 /* CRBIV1: Table of coeffs of the approximation of iso-V0 and its */
1880 /* derivatives till order IORDRV. */
1881 /* NCFIV2: Nb of Coeff. of curves stored in CRBIV2. */
1882 /* CRBIV2: Table of coeffs of approximation of iso-V1 and its */
1883 /* derivatives till order IORDRV. */
1884 /* VHERMT: Coeff. of Hermit polynoms of order IORDRV. */
1885 /* PATJAC: Table of coefficients of the polynom P(u,v) of approximation */
1886 /* of F(u,v) WITHOUT taking into account the constraints. */
1888 /* OUTPUT ARGUMENTS : */
1889 /* ------------------- */
1890 /* PATJAC: Table of coefficients of the polynom P(u,v) by approximation */
1891 /* of F(u,v) WITH taking into account of constraints. */
1896 /* **********************************************************************
1898 /* Name of the routine */
1900 /* --------------------------- Initialisations --------------------------
1903 /* Parameter adjustments */
1904 patjac_dim1 = *mxujac + 1;
1905 patjac_dim2 = *mxvjac + 1;
1906 patjac_offset = patjac_dim1 * patjac_dim2;
1907 patjac -= patjac_offset;
1908 vhermt_dim1 = (*iordrv << 1) + 2;
1909 vhermt_offset = vhermt_dim1;
1910 vhermt -= vhermt_offset;
1913 crbiv2_dim1 = *nclimu;
1914 crbiv2_dim2 = *ndimen;
1915 crbiv2_offset = crbiv2_dim1 * (crbiv2_dim2 + 1);
1916 crbiv2 -= crbiv2_offset;
1917 crbiv1_dim1 = *nclimu;
1918 crbiv1_dim2 = *ndimen;
1919 crbiv1_offset = crbiv1_dim1 * (crbiv1_dim2 + 1);
1920 crbiv1 -= crbiv1_offset;
1923 ldbg = AdvApp2Var_SysBase::mnfndeb_() >= 3;
1925 AdvApp2Var_SysBase::mgenmsg_("MMA2AC2", 7L);
1928 /* ------------ ADDING of coeff by u of curves, by v of Hermit --------
1932 for (ii = 1; ii <= i__1; ++ii) {
1933 ndgv1 = ncfiv1[ii] - 1;
1934 ndgv2 = ncfiv2[ii] - 1;
1936 for (nd = 1; nd <= i__2; ++nd) {
1937 i__3 = (*iordrv << 1) + 1;
1938 for (jj = 0; jj <= i__3; ++jj) {
1939 bid1 = vhermt[jj + ((ii << 1) - 1) * vhermt_dim1];
1941 for (kk = 0; kk <= i__4; ++kk) {
1942 patjac[kk + (jj + nd * patjac_dim2) * patjac_dim1] +=
1943 bid1 * crbiv1[kk + (nd + ii * crbiv1_dim2) *
1947 bid2 = vhermt[jj + (ii << 1) * vhermt_dim1];
1949 for (kk = 0; kk <= i__4; ++kk) {
1950 patjac[kk + (jj + nd * patjac_dim2) * patjac_dim1] +=
1951 bid2 * crbiv2[kk + (nd + ii * crbiv2_dim2) *
1962 /* ------------------------------ The end -------------------------------
1966 AdvApp2Var_SysBase::mgsomsg_("MMA2AC2", 7L);
1972 //=======================================================================
1973 //function : mma2ac3_
1975 //=======================================================================
1976 int AdvApp2Var_ApproxF2var::mma2ac3_(const integer *ndimen,
1977 const integer *mxujac,
1978 const integer *mxvjac,
1979 const integer *iordru,
1980 const integer *nclimv,
1981 const integer *ncfiu1,
1982 const doublereal * crbiu1,
1983 const integer *ncfiu2,
1984 const doublereal *crbiu2,
1985 const doublereal *uhermt,
1989 /* System generated locals */
1990 integer crbiu1_dim1, crbiu1_dim2, crbiu1_offset, crbiu2_dim1, crbiu2_dim2,
1991 crbiu2_offset, patjac_dim1, patjac_dim2, patjac_offset,
1992 uhermt_dim1, uhermt_offset, i__1, i__2, i__3, i__4;
1994 /* Local variables */
1995 static logical ldbg;
1996 static integer ndgu1, ndgu2, ii, jj, nd, kk;
1997 static doublereal bid1, bid2;
1999 /* **********************************************************************
2004 /* Ajout des polynomes de contraintes */
2008 /* FONCTION,APPROXIMATION,COEFFICIENT,POLYNOME */
2010 /* INPUT ARGUMENTS : */
2011 /* ------------------ */
2012 /* NDIMEN: Dimension of the space. */
2013 /* MXUJAC: Max degree of the polynom of approximation by U. The */
2014 /* representation in the orthogonal base starts from degree */
2015 /* 0 to degree MXUJAC-2*(IORDRU+1). The polynomial base is the */
2016 /* base of Jacobi of order -1 (Legendre), 0, 1 or 2 */
2017 /* MXVJAC: Max degree of the polynom of approximation by V. The */
2018 /* representation in the orthogonal base starts from degree */
2019 /* 0 to degree MXUJAC-2*(IORDRU+1). The polynomial base is the */
2020 /* base of Jacobi of order -1 (Legendre), 0, 1 or 2 */
2021 /* IORDRU: Order of the base of Jacobi (-1,0,1 or 2) by U. Corresponds */
2022 /* to the step of constraints: C0, C1 or C2. */
2023 /* NCLIMV LIMIT nb of coeff by v of the solution P(u,v)
2024 * NCFIU1: Nb of Coeff. of curves stored in CRBIU1. */
2025 /* CRBIU1: Table of coeffs of the approximation of iso-U0 and its */
2026 /* derivatives till order IORDRU. */
2027 /* NCFIU2: Nb of Coeff. of curves stored in CRBIU2. */
2028 /* CRBIU2: Table of coeffs of approximation of iso-U1 and its */
2029 /* derivatives till order IORDRU */
2030 /* UHERMT: Coeff. of Hermit polynoms of order IORDRU. */
2031 /* PATJAC: Table of coefficients of the polynom P(u,v) of approximation */
2032 /* of F(u,v) WITHOUT taking into account the constraints. */
2034 /* OUTPUT ARGUMENTS : */
2035 /* ------------------- */
2036 /* PATJAC: Table of coefficients of the polynom P(u,v) by approximation */
2037 /* of F(u,v) WITH taking into account of constraints. */
2041 /* **********************************************************************
2043 /* The name of the routine */
2045 /* --------------------------- Initializations --------------------------
2048 /* Parameter adjustments */
2049 patjac_dim1 = *mxujac + 1;
2050 patjac_dim2 = *mxvjac + 1;
2051 patjac_offset = patjac_dim1 * patjac_dim2;
2052 patjac -= patjac_offset;
2053 uhermt_dim1 = (*iordru << 1) + 2;
2054 uhermt_offset = uhermt_dim1;
2055 uhermt -= uhermt_offset;
2058 crbiu2_dim1 = *nclimv;
2059 crbiu2_dim2 = *ndimen;
2060 crbiu2_offset = crbiu2_dim1 * (crbiu2_dim2 + 1);
2061 crbiu2 -= crbiu2_offset;
2062 crbiu1_dim1 = *nclimv;
2063 crbiu1_dim2 = *ndimen;
2064 crbiu1_offset = crbiu1_dim1 * (crbiu1_dim2 + 1);
2065 crbiu1 -= crbiu1_offset;
2068 ldbg = AdvApp2Var_SysBase::mnfndeb_() >= 3;
2070 AdvApp2Var_SysBase::mgenmsg_("MMA2AC3", 7L);
2073 /* ------------ ADDING of coeff by u of curves, by v of Hermit --------
2077 for (ii = 1; ii <= i__1; ++ii) {
2078 ndgu1 = ncfiu1[ii] - 1;
2079 ndgu2 = ncfiu2[ii] - 1;
2081 for (nd = 1; nd <= i__2; ++nd) {
2083 for (jj = 0; jj <= i__3; ++jj) {
2084 bid1 = crbiu1[jj + (nd + ii * crbiu1_dim2) * crbiu1_dim1];
2085 i__4 = (*iordru << 1) + 1;
2086 for (kk = 0; kk <= i__4; ++kk) {
2087 patjac[kk + (jj + nd * patjac_dim2) * patjac_dim1] +=
2088 bid1 * uhermt[kk + ((ii << 1) - 1) * uhermt_dim1];
2094 for (jj = 0; jj <= i__3; ++jj) {
2095 bid2 = crbiu2[jj + (nd + ii * crbiu2_dim2) * crbiu2_dim1];
2096 i__4 = (*iordru << 1) + 1;
2097 for (kk = 0; kk <= i__4; ++kk) {
2098 patjac[kk + (jj + nd * patjac_dim2) * patjac_dim1] +=
2099 bid2 * uhermt[kk + (ii << 1) * uhermt_dim1];
2110 /* ------------------------------ The end -------------------------------
2114 AdvApp2Var_SysBase::mgsomsg_("MMA2AC3", 7L);
2119 //=======================================================================
2120 //function : mma2can_
2122 //=======================================================================
2123 int AdvApp2Var_ApproxF2var::mma2can_(const integer *ncfmxu,
2124 const integer *ncfmxv,
2125 const integer *ndimen,
2126 const integer *iordru,
2127 const integer *iordrv,
2128 const integer *ncoefu,
2129 const integer *ncoefv,
2130 const doublereal *patjac,
2136 /* System generated locals */
2137 integer patjac_dim1, patjac_dim2, patjac_offset, patcan_dim1, patcan_dim2,
2138 patcan_offset, i__1, i__2;
2140 /* Local variables */
2141 static logical ldbg;
2142 static integer ilon1, ilon2, ii, nd;
2144 /* **********************************************************************
2149 /* Change of Jacobi base to canonical (-1,1) and writing in a greater */
2154 /* ALL,AB_SPECIFI,CARREAU&,CONVERSION,JACOBI,CANNONIQUE,&CARREAU */
2156 /* INPUT ARGUMENTS : */
2157 /* -------------------- */
2158 /* NCFMXU: Dimension by U of resulting table PATCAN */
2159 /* NCFMXV: Dimension by V of resulting table PATCAN */
2160 /* NDIMEN: Dimension of the workspace. */
2161 /* IORDRU: Order of constraint by U */
2162 /* IORDRV: Order of constraint by V. */
2163 /* NCOEFU: Nb of coeff by U of square PATJAC */
2164 /* NCOEFV: Nb of coeff by V of square PATJAC */
2165 /* PATJAC: Square in the base of Jacobi of order IORDRU by U and */
2168 /* OUTPUT ARGUMENTS : */
2169 /* --------------------- */
2170 /* PATAUX: Auxiliary Table. */
2171 /* PATCAN: Table of coefficients in the canonic base. */
2172 /* IERCOD: Error code. */
2173 /* = 0, everything goes well, and all things are equal. */
2174 /* = 1, the program refuses to process with incorrect input arguments */
2177 /* COMMONS USED : */
2178 /* ------------------ */
2180 /* REFERENCES CALLED : */
2181 /* --------------------- */
2183 /* DESCRIPTION/NOTES/LIMITATIONS : */
2184 /* ----------------------------------- */
2186 /* **********************************************************************
2190 /* Parameter adjustments */
2191 patcan_dim1 = *ncfmxu;
2192 patcan_dim2 = *ncfmxv;
2193 patcan_offset = patcan_dim1 * (patcan_dim2 + 1) + 1;
2194 patcan -= patcan_offset;
2196 patjac_dim1 = *ncoefu;
2197 patjac_dim2 = *ncoefv;
2198 patjac_offset = patjac_dim1 * (patjac_dim2 + 1) + 1;
2199 patjac -= patjac_offset;
2202 ldbg = AdvApp2Var_SysBase::mnfndeb_() >= 2;
2204 AdvApp2Var_SysBase::mgenmsg_("MMA2CAN", 7L);
2208 if (*iordru < -1 || *iordru > 2) {
2211 if (*iordrv < -1 || *iordrv > 2) {
2214 if (*ncoefu > *ncfmxu || *ncoefv > *ncfmxv) {
2218 /* --> Pass to canonic base (-1,1) */
2219 mmjacpt_(ndimen, ncoefu, ncoefv, iordru, iordrv, &patjac[patjac_offset], &
2220 pataux[1], &patcan[patcan_offset]);
2222 /* --> Write all in a greater table */
2223 AdvApp2Var_MathBase::mmfmca8_((integer *)ncoefu,
2229 (doublereal *)&patcan[patcan_offset],
2230 (doublereal *)&patcan[patcan_offset]);
2232 /* --> Complete with zeros the resulting table. */
2233 ilon1 = *ncfmxu - *ncoefu;
2234 ilon2 = *ncfmxu * (*ncfmxv - *ncoefv);
2236 for (nd = 1; nd <= i__1; ++nd) {
2239 for (ii = 1; ii <= i__2; ++ii) {
2240 AdvApp2Var_SysBase::mvriraz_(&ilon1,
2241 (char *)&patcan[*ncoefu + 1 + (ii + nd * patcan_dim2) * patcan_dim1]);
2246 AdvApp2Var_SysBase::mvriraz_(&ilon2,
2247 (char *)&patcan[(*ncoefv + 1 + nd * patcan_dim2) * patcan_dim1 + 1]);
2254 /* ----------------------
2262 AdvApp2Var_SysBase::maermsg_("MMA2CAN", iercod, 7L);
2264 AdvApp2Var_SysBase::mgsomsg_("MMA2CAN", 7L);
2269 //=======================================================================
2270 //function : mma2cd1_
2272 //=======================================================================
2273 int mma2cd1_(integer *ndimen,
2294 static integer c__1 = 1;
2296 /* System generated locals */
2297 integer contr1_dim1, contr1_dim2, contr1_offset, contr2_dim1, contr2_dim2,
2298 contr2_offset, contr3_dim1, contr3_dim2, contr3_offset,
2299 contr4_dim1, contr4_dim2, contr4_offset, uhermt_dim1,
2300 uhermt_offset, vhermt_dim1, vhermt_offset, fpntbu_dim1,
2301 fpntbu_offset, fpntbv_dim1, fpntbv_offset, sosotb_dim1,
2302 sosotb_dim2, sosotb_offset, diditb_dim1, diditb_dim2,
2303 diditb_offset, soditb_dim1, soditb_dim2, soditb_offset,
2304 disotb_dim1, disotb_dim2, disotb_offset, i__1, i__2, i__3, i__4,
2307 /* Local variables */
2308 static integer ncfhu, ncfhv, nuroo, nvroo, nd, ii, jj, kk, ll, ibb, kkm,
2310 static doublereal bid1, bid2, bid3, bid4;
2311 static doublereal diu1, diu2, div1, div2, sou1, sou2, sov1, sov2;
2313 /* **********************************************************************
2318 /* Discretisation on the parameters of polynoms of interpolation */
2319 /* of constraints at the corners of order IORDRE. */
2323 /* TOUS, AB_SPECIFI::CONTRAINTE&, DISCRETISATION, &POINT */
2325 /* INPUT ARGUMENTS : */
2326 /* ------------------ */
2327 /* NDIMEN: Dimension of the space. */
2328 /* NBPNTU: Nb of INTERNAL parameters of discretisation by U. */
2329 /* This is also the nb of root of Legendre polynom where discretization is done. */
2330 /* UROOTL: Table of parameters of discretisation ON (-1,1) by U.
2332 /* NBPNTV: Nb of INTERNAL parameters of discretisation by V. */
2333 /* This is also the nb of root of Legendre polynom where discretization is done. */
2334 /* VROOTL: Table of discretization parameters on (-1,1) by V.
2335 /* IORDRU: Order of constraint imposed at the extremities of iso-V */
2336 /* = 0, calculate the extremities of iso-V */
2337 /* = 1, calculate, additionally, the 1st derivative in the direction of iso-V */
2338 /* = 2, calculate, additionally, the 2nd derivative in the direction of iso-V */
2339 /* IORDRV: Order of constraint imposed at the extremities of iso-U */
2340 /* = 0, calculate the extremities of iso-U */
2341 /* = 1, calculate, additionally, the 1st derivative in the direction of iso-U */
2342 /* = 2, calculate, additionally, the 2nd derivative in the direction of iso-U */
2343 /* CONTR1: Contains, if IORDRU and IORDRV>=0, the values at the */
2344 /* extremities of F(U0,V0) and its derivatives. */
2345 /* CONTR2: Contains, if IORDRU and IORDRV>=0, the values at the */
2346 /* extremities of F(U1,V0) and its derivatives. */
2347 /* CONTR3: Contains, if IORDRU and IORDRV>=0, the values at the */
2348 /* extremities of F(U0,V1) and its derivatives. */
2349 /* CONTR4: Contains, if IORDRU and IORDRV>=0, the values at the */
2350 /* extremities of F(U1,V1) and its derivatives. */
2351 /* SOSOTB: Preinitialized table (input/output argument). */
2352 /* DISOTB: Preinitialized table (input/output argument). */
2353 /* SODITB: Preinitialized table (input/output argument). */
2354 /* DIDITB: Preinitialized table (input/output argument) */
2356 /* OUTPUT ARGUMENTS : */
2357 /* ------------------- */
2358 /* FPNTBU: Auxiliary table. */
2359 /* FPNTBV: Auxiliary table. */
2360 /* UHERMT: Table of 2*(IORDRU+1) coeff. of 2*(IORDRU+1) polynoms of Hermite. */
2361 /* VHERMT: Table of 2*(IORDRV+1) coeff. of 2*(IORDRV+1) polynoms of Hermite. */
2362 /* SOSOTB: Table where the terms of constraints are added */
2363 /* C(ui,vj) + C(ui,-vj) + C(-ui,vj) + C(-ui,-vj) */
2364 /* with ui and vj positive roots of the Legendre polynom */
2365 /* of degree NBPNTU and NBPNTV respectively. */
2366 /* DISOTB: Table where the terms of constraints are added */
2367 /* C(ui,vj) + C(ui,-vj) - C(-ui,vj) - C(-ui,-vj) */
2368 /* with ui and vj positive roots of the polynom of Legendre */
2369 /* of degree NBPNTU and NBPNTV respectively. */
2370 /* SODITB: Table where the terms of constraints are added */
2371 /* C(ui,vj) - C(ui,-vj) + C(-ui,vj) - C(-ui,-vj) */
2372 /* with ui and vj positive roots of the polynom of Legendre */
2373 /* of degree NBPNTU and NBPNTV respectively. */
2374 /* DIDITB: Table where the terms of constraints are added */
2375 /* C(ui,vj) - C(ui,-vj) - C(-ui,vj) + C(-ui,-vj) */
2376 /* with ui and vj positive roots of the polynom of Legendre */
2377 /* of degree NBPNTU and NBPNTV respectively. */
2379 /* COMMONS USED : */
2380 /* ---------------- */
2382 /* REFERENCES CALLED : */
2383 /* ----------------------- */
2385 /* DESCRIPTION/NOTES/LIMITATIONS : */
2386 /* ----------------------------------- */
2389 /* **********************************************************************
2392 /* Name of the routine */
2395 /* Parameter adjustments */
2397 diditb_dim1 = *nbpntu / 2 + 1;
2398 diditb_dim2 = *nbpntv / 2 + 1;
2399 diditb_offset = diditb_dim1 * diditb_dim2;
2400 diditb -= diditb_offset;
2401 disotb_dim1 = *nbpntu / 2;
2402 disotb_dim2 = *nbpntv / 2;
2403 disotb_offset = disotb_dim1 * (disotb_dim2 + 1) + 1;
2404 disotb -= disotb_offset;
2405 soditb_dim1 = *nbpntu / 2;
2406 soditb_dim2 = *nbpntv / 2;
2407 soditb_offset = soditb_dim1 * (soditb_dim2 + 1) + 1;
2408 soditb -= soditb_offset;
2409 sosotb_dim1 = *nbpntu / 2 + 1;
2410 sosotb_dim2 = *nbpntv / 2 + 1;
2411 sosotb_offset = sosotb_dim1 * sosotb_dim2;
2412 sosotb -= sosotb_offset;
2414 uhermt_dim1 = (*iordru << 1) + 2;
2415 uhermt_offset = uhermt_dim1;
2416 uhermt -= uhermt_offset;
2417 fpntbu_dim1 = *nbpntu;
2418 fpntbu_offset = fpntbu_dim1 + 1;
2419 fpntbu -= fpntbu_offset;
2420 vhermt_dim1 = (*iordrv << 1) + 2;
2421 vhermt_offset = vhermt_dim1;
2422 vhermt -= vhermt_offset;
2423 fpntbv_dim1 = *nbpntv;
2424 fpntbv_offset = fpntbv_dim1 + 1;
2425 fpntbv -= fpntbv_offset;
2426 contr4_dim1 = *ndimen;
2427 contr4_dim2 = *iordru + 2;
2428 contr4_offset = contr4_dim1 * (contr4_dim2 + 1) + 1;
2429 contr4 -= contr4_offset;
2430 contr3_dim1 = *ndimen;
2431 contr3_dim2 = *iordru + 2;
2432 contr3_offset = contr3_dim1 * (contr3_dim2 + 1) + 1;
2433 contr3 -= contr3_offset;
2434 contr2_dim1 = *ndimen;
2435 contr2_dim2 = *iordru + 2;
2436 contr2_offset = contr2_dim1 * (contr2_dim2 + 1) + 1;
2437 contr2 -= contr2_offset;
2438 contr1_dim1 = *ndimen;
2439 contr1_dim2 = *iordru + 2;
2440 contr1_offset = contr1_dim1 * (contr1_dim2 + 1) + 1;
2441 contr1 -= contr1_offset;
2444 ibb = AdvApp2Var_SysBase::mnfndeb_();
2446 AdvApp2Var_SysBase::mgenmsg_("MMA2CD1", 7L);
2449 /* ------------------- Discretisation of Hermite polynoms -----------
2452 ncfhu = (*iordru + 1) << 1;
2454 for (ii = 1; ii <= i__1; ++ii) {
2456 for (ll = 1; ll <= i__2; ++ll) {
2457 AdvApp2Var_MathBase::mmmpocur_(&ncfhu, &c__1, &ncfhu, &uhermt[ii * uhermt_dim1], &
2458 urootl[ll], &fpntbu[ll + ii * fpntbu_dim1]);
2463 ncfhv = (*iordrv + 1) << 1;
2465 for (jj = 1; jj <= i__1; ++jj) {
2467 for (kk = 1; kk <= i__2; ++kk) {
2468 AdvApp2Var_MathBase::mmmpocur_(&ncfhv, &c__1, &ncfhv, &vhermt[jj * vhermt_dim1], &
2469 vrootl[kk], &fpntbv[kk + jj * fpntbv_dim1]);
2475 /* ---- The discretizations of polynoms of constraints are subtracted ----
2478 nuroo = *nbpntu / 2;
2479 nvroo = *nbpntv / 2;
2481 for (nd = 1; nd <= i__1; ++nd) {
2484 for (jj = 1; jj <= i__2; ++jj) {
2486 for (ii = 1; ii <= i__3; ++ii) {
2487 bid1 = contr1[nd + (ii + jj * contr1_dim2) * contr1_dim1];
2488 bid2 = contr2[nd + (ii + jj * contr2_dim2) * contr2_dim1];
2489 bid3 = contr3[nd + (ii + jj * contr3_dim2) * contr3_dim1];
2490 bid4 = contr4[nd + (ii + jj * contr4_dim2) * contr4_dim1];
2493 for (kk = 1; kk <= i__4; ++kk) {
2494 kkp = (*nbpntv + 1) / 2 + kk;
2495 kkm = nvroo - kk + 1;
2496 sov1 = fpntbv[kkp + ((jj << 1) - 1) * fpntbv_dim1] +
2497 fpntbv[kkm + ((jj << 1) - 1) * fpntbv_dim1];
2498 div1 = fpntbv[kkp + ((jj << 1) - 1) * fpntbv_dim1] -
2499 fpntbv[kkm + ((jj << 1) - 1) * fpntbv_dim1];
2500 sov2 = fpntbv[kkp + (jj << 1) * fpntbv_dim1] + fpntbv[kkm
2501 + (jj << 1) * fpntbv_dim1];
2502 div2 = fpntbv[kkp + (jj << 1) * fpntbv_dim1] - fpntbv[kkm
2503 + (jj << 1) * fpntbv_dim1];
2505 for (ll = 1; ll <= i__5; ++ll) {
2506 llp = (*nbpntu + 1) / 2 + ll;
2507 llm = nuroo - ll + 1;
2508 sou1 = fpntbu[llp + ((ii << 1) - 1) * fpntbu_dim1] +
2509 fpntbu[llm + ((ii << 1) - 1) * fpntbu_dim1];
2510 diu1 = fpntbu[llp + ((ii << 1) - 1) * fpntbu_dim1] -
2511 fpntbu[llm + ((ii << 1) - 1) * fpntbu_dim1];
2512 sou2 = fpntbu[llp + (ii << 1) * fpntbu_dim1] + fpntbu[
2513 llm + (ii << 1) * fpntbu_dim1];
2514 diu2 = fpntbu[llp + (ii << 1) * fpntbu_dim1] - fpntbu[
2515 llm + (ii << 1) * fpntbu_dim1];
2516 sosotb[ll + (kk + nd * sosotb_dim2) * sosotb_dim1] =
2517 sosotb[ll + (kk + nd * sosotb_dim2) *
2518 sosotb_dim1] - bid1 * sou1 * sov1 - bid2 *
2519 sou2 * sov1 - bid3 * sou1 * sov2 - bid4 *
2521 soditb[ll + (kk + nd * soditb_dim2) * soditb_dim1] =
2522 soditb[ll + (kk + nd * soditb_dim2) *
2523 soditb_dim1] - bid1 * sou1 * div1 - bid2 *
2524 sou2 * div1 - bid3 * sou1 * div2 - bid4 *
2526 disotb[ll + (kk + nd * disotb_dim2) * disotb_dim1] =
2527 disotb[ll + (kk + nd * disotb_dim2) *
2528 disotb_dim1] - bid1 * diu1 * sov1 - bid2 *
2529 diu2 * sov1 - bid3 * diu1 * sov2 - bid4 *
2531 diditb[ll + (kk + nd * diditb_dim2) * diditb_dim1] =
2532 diditb[ll + (kk + nd * diditb_dim2) *
2533 diditb_dim1] - bid1 * diu1 * div1 - bid2 *
2534 diu2 * div1 - bid3 * diu1 * div2 - bid4 *
2541 /* ------------ Case when the discretization is done only on the roots
2543 /* ---------- of Legendre polynom of uneven degree, 0 is root
2546 if (*nbpntu % 2 == 1) {
2547 sou1 = fpntbu[nuroo + 1 + ((ii << 1) - 1) * fpntbu_dim1];
2548 sou2 = fpntbu[nuroo + 1 + (ii << 1) * fpntbu_dim1];
2550 for (kk = 1; kk <= i__4; ++kk) {
2551 kkp = (*nbpntv + 1) / 2 + kk;
2552 kkm = nvroo - kk + 1;
2553 sov1 = fpntbv[kkp + ((jj << 1) - 1) * fpntbv_dim1] +
2554 fpntbv[kkm + ((jj << 1) - 1) * fpntbv_dim1];
2555 div1 = fpntbv[kkp + ((jj << 1) - 1) * fpntbv_dim1] -
2556 fpntbv[kkm + ((jj << 1) - 1) * fpntbv_dim1];
2557 sov2 = fpntbv[kkp + (jj << 1) * fpntbv_dim1] + fpntbv[
2558 kkm + (jj << 1) * fpntbv_dim1];
2559 div2 = fpntbv[kkp + (jj << 1) * fpntbv_dim1] - fpntbv[
2560 kkm + (jj << 1) * fpntbv_dim1];
2561 sosotb[(kk + nd * sosotb_dim2) * sosotb_dim1] =
2562 sosotb[(kk + nd * sosotb_dim2) * sosotb_dim1]
2563 - bid1 * sou1 * sov1 - bid2 * sou2 * sov1 -
2564 bid3 * sou1 * sov2 - bid4 * sou2 * sov2;
2565 diditb[(kk + nd * diditb_dim2) * diditb_dim1] =
2566 diditb[(kk + nd * diditb_dim2) * diditb_dim1]
2567 - bid1 * sou1 * div1 - bid2 * sou2 * div1 -
2568 bid3 * sou1 * div2 - bid4 * sou2 * div2;
2573 if (*nbpntv % 2 == 1) {
2574 sov1 = fpntbv[nvroo + 1 + ((jj << 1) - 1) * fpntbv_dim1];
2575 sov2 = fpntbv[nvroo + 1 + (jj << 1) * fpntbv_dim1];
2577 for (ll = 1; ll <= i__4; ++ll) {
2578 llp = (*nbpntu + 1) / 2 + ll;
2579 llm = nuroo - ll + 1;
2580 sou1 = fpntbu[llp + ((ii << 1) - 1) * fpntbu_dim1] +
2581 fpntbu[llm + ((ii << 1) - 1) * fpntbu_dim1];
2582 diu1 = fpntbu[llp + ((ii << 1) - 1) * fpntbu_dim1] -
2583 fpntbu[llm + ((ii << 1) - 1) * fpntbu_dim1];
2584 sou2 = fpntbu[llp + (ii << 1) * fpntbu_dim1] + fpntbu[
2585 llm + (ii << 1) * fpntbu_dim1];
2586 diu2 = fpntbu[llp + (ii << 1) * fpntbu_dim1] - fpntbu[
2587 llm + (ii << 1) * fpntbu_dim1];
2588 sosotb[ll + nd * sosotb_dim2 * sosotb_dim1] = sosotb[
2589 ll + nd * sosotb_dim2 * sosotb_dim1] - bid1 *
2590 sou1 * sov1 - bid2 * sou2 * sov1 - bid3 *
2591 sou1 * sov2 - bid4 * sou2 * sov2;
2592 diditb[ll + nd * diditb_dim2 * diditb_dim1] = diditb[
2593 ll + nd * diditb_dim2 * diditb_dim1] - bid1 *
2594 diu1 * sov1 - bid2 * diu2 * sov1 - bid3 *
2595 diu1 * sov2 - bid4 * diu2 * sov2;
2600 if (*nbpntu % 2 == 1 && *nbpntv % 2 == 1) {
2601 sou1 = fpntbu[nuroo + 1 + ((ii << 1) - 1) * fpntbu_dim1];
2602 sou2 = fpntbu[nuroo + 1 + (ii << 1) * fpntbu_dim1];
2603 sov1 = fpntbv[nvroo + 1 + ((jj << 1) - 1) * fpntbv_dim1];
2604 sov2 = fpntbv[nvroo + 1 + (jj << 1) * fpntbv_dim1];
2605 sosotb[nd * sosotb_dim2 * sosotb_dim1] = sosotb[nd *
2606 sosotb_dim2 * sosotb_dim1] - bid1 * sou1 * sov1 -
2607 bid2 * sou2 * sov1 - bid3 * sou1 * sov2 - bid4 *
2609 diditb[nd * diditb_dim2 * diditb_dim1] = diditb[nd *
2610 diditb_dim2 * diditb_dim1] - bid1 * sou1 * sov1 -
2611 bid2 * sou2 * sov1 - bid3 * sou1 * sov2 - bid4 *
2623 /* ------------------------------ The End -------------------------------
2628 AdvApp2Var_SysBase::mgsomsg_("MMA2CD1", 7L);
2633 //=======================================================================
2634 //function : mma2cd2_
2636 //=======================================================================
2637 int mma2cd2_(integer *ndimen,
2654 static integer c__1 = 1;
2655 /* System generated locals */
2656 integer sotbv1_dim1, sotbv1_dim2, sotbv1_offset, sotbv2_dim1, sotbv2_dim2,
2657 sotbv2_offset, ditbv1_dim1, ditbv1_dim2, ditbv1_offset,
2658 ditbv2_dim1, ditbv2_dim2, ditbv2_offset, fpntab_dim1,
2659 fpntab_offset, vhermt_dim1, vhermt_offset, sosotb_dim1,
2660 sosotb_dim2, sosotb_offset, diditb_dim1, diditb_dim2,
2661 diditb_offset, soditb_dim1, soditb_dim2, soditb_offset,
2662 disotb_dim1, disotb_dim2, disotb_offset, i__1, i__2, i__3, i__4;
2664 /* Local variables */
2665 static integer ncfhv, nuroo, nvroo, ii, nd, jj, kk, ibb, jjm, jjp;
2666 static doublereal bid1, bid2, bid3, bid4;
2668 /* **********************************************************************
2672 /* Discretisation on the parameters of polynoms of interpolation */
2673 /* of constraints on 2 borders iso-V of order IORDRV. */
2678 /* TOUS, AB_SPECIFI::CONTRAINTE&, DISCRETISATION, &POINT */
2682 /* INPUT ARGUMENTS : */
2683 /* ------------------ */
2684 /* NDIMEN: Dimension of the space. */
2685 /* NBPNTU: Nb of INTERNAL parameters of discretisation by U. */
2686 /* This is also the nb of root of Legendre polynom where discretization is done. */
2687 /* UROOTL: Table of parameters of discretisation ON (-1,1) by U.
2689 /* NBPNTV: Nb of INTERNAL parameters of discretisation by V. */
2690 /* This is also the nb of root of Legendre polynom where discretization is done. */
2691 /* VROOTL: Table of discretization parameters on (-1,1) by V.
2692 /* IORDRV: Order of constraint imposed at the extremities of iso-V */
2693 /* = 0, calculate the extremities of iso-V */
2694 /* = 1, calculate, additionally, the 1st derivative in the direction of iso-V */
2695 /* = 2, calculate, additionally, the 2nd derivative in the direction of iso-V */
2696 /* SOTBV1: Table of NBPNTV/2 sums of 2 index points */
2697 /* NBPNTV-II+1 and II, for II = 1, NBPNTV/2 on iso-V0. */
2698 /* SOTBV2: Table of NBPNTV/2 sums of 2 index points */
2699 /* NBPNTV-II+1 and II, for II = 1, NBPNTV/2 on iso-V1. */
2700 /* DITBV1: Table of NBPNTV/2 differences of 2 index points */
2701 /* NBPNTV-II+1 and II, for II = 1, NBPNTV/2 on iso-V0. */
2702 /* DITBV2: Table of NBPNTV/2 differences of 2 index points */
2703 /* NBPNTV-II+1 and II, for II = 1, NBPNTV/2 on iso-V1. */
2704 /* SOSOTB: Preinitialized table (input/output argument). */
2705 /* DISOTB: Preinitialized table (input/output argument). */
2706 /* SODITB: Preinitialized table (input/output argument). */
2707 /* DIDITB: Preinitialized table (input/output argument) */
2709 /* OUTPUT ARGUMENTS : */
2710 /* ------------------- */
2711 /* FPNTAB: Auxiliary table. */
2712 /* VHERMT: Table of 2*(IORDRV+1) coeff. of 2*(IORDRV+1) polynoms of Hermite. */
2713 /* SOSOTB: Table where the terms of constraints are added */
2714 /* C(ui,vj) + C(ui,-vj) + C(-ui,vj) + C(-ui,-vj) */
2715 /* with ui and vj positive roots of the Legendre polynom */
2716 /* of degree NBPNTU and NBPNTV respectively. */
2717 /* DISOTB: Table where the terms of constraints are added */
2718 /* C(ui,vj) + C(ui,-vj) - C(-ui,vj) - C(-ui,-vj) */
2719 /* with ui and vj positive roots of the polynom of Legendre */
2720 /* of degree NBPNTU and NBPNTV respectively. */
2721 /* SODITB: Table where the terms of constraints are added */
2722 /* C(ui,vj) - C(ui,-vj) + C(-ui,vj) - C(-ui,-vj) */
2723 /* with ui and vj positive roots of the polynom of Legendre */
2724 /* of degree NBPNTU and NBPNTV respectively. */
2725 /* DIDITB: Table where the terms of constraints are added */
2726 /* C(ui,vj) - C(ui,-vj) - C(-ui,vj) + C(-ui,-vj) */
2727 /* with ui and vj positive roots of the polynom of Legendre */
2728 /* of degree NBPNTU and NBPNTV respectively. */
2730 /* COMMONS USED : */
2731 /* ---------------- */
2733 /* REFERENCES CALLED : */
2734 /* ----------------------- */
2736 /* DESCRIPTION/NOTES/LIMITATIONS : */
2737 /* ----------------------------------- */
2741 /* **********************************************************************
2744 /* Name of the routine */
2747 /* Parameter adjustments */
2748 diditb_dim1 = *nbpntu / 2 + 1;
2749 diditb_dim2 = *nbpntv / 2 + 1;
2750 diditb_offset = diditb_dim1 * diditb_dim2;
2751 diditb -= diditb_offset;
2752 disotb_dim1 = *nbpntu / 2;
2753 disotb_dim2 = *nbpntv / 2;
2754 disotb_offset = disotb_dim1 * (disotb_dim2 + 1) + 1;
2755 disotb -= disotb_offset;
2756 soditb_dim1 = *nbpntu / 2;
2757 soditb_dim2 = *nbpntv / 2;
2758 soditb_offset = soditb_dim1 * (soditb_dim2 + 1) + 1;
2759 soditb -= soditb_offset;
2760 sosotb_dim1 = *nbpntu / 2 + 1;
2761 sosotb_dim2 = *nbpntv / 2 + 1;
2762 sosotb_offset = sosotb_dim1 * sosotb_dim2;
2763 sosotb -= sosotb_offset;
2765 vhermt_dim1 = (*iordrv << 1) + 2;
2766 vhermt_offset = vhermt_dim1;
2767 vhermt -= vhermt_offset;
2768 fpntab_dim1 = *nbpntv;
2769 fpntab_offset = fpntab_dim1 + 1;
2770 fpntab -= fpntab_offset;
2771 ditbv2_dim1 = *nbpntu / 2 + 1;
2772 ditbv2_dim2 = *ndimen;
2773 ditbv2_offset = ditbv2_dim1 * (ditbv2_dim2 + 1);
2774 ditbv2 -= ditbv2_offset;
2775 ditbv1_dim1 = *nbpntu / 2 + 1;
2776 ditbv1_dim2 = *ndimen;
2777 ditbv1_offset = ditbv1_dim1 * (ditbv1_dim2 + 1);
2778 ditbv1 -= ditbv1_offset;
2779 sotbv2_dim1 = *nbpntu / 2 + 1;
2780 sotbv2_dim2 = *ndimen;
2781 sotbv2_offset = sotbv2_dim1 * (sotbv2_dim2 + 1);
2782 sotbv2 -= sotbv2_offset;
2783 sotbv1_dim1 = *nbpntu / 2 + 1;
2784 sotbv1_dim2 = *ndimen;
2785 sotbv1_offset = sotbv1_dim1 * (sotbv1_dim2 + 1);
2786 sotbv1 -= sotbv1_offset;
2789 ibb = AdvApp2Var_SysBase::mnfndeb_();
2791 AdvApp2Var_SysBase::mgenmsg_("MMA2CD2", 7L);
2794 /* ------------------- Discretization of Hermit polynoms -----------
2797 ncfhv = (*iordrv + 1) << 1;
2799 for (ii = 1; ii <= i__1; ++ii) {
2801 for (jj = 1; jj <= i__2; ++jj) {
2802 AdvApp2Var_MathBase::mmmpocur_(&ncfhv, &c__1, &ncfhv, &vhermt[ii * vhermt_dim1], &
2803 vrootl[jj], &fpntab[jj + ii * fpntab_dim1]);
2809 /* ---- The discretizations of polynoms of constraints are subtracted ----
2812 nuroo = *nbpntu / 2;
2813 nvroo = *nbpntv / 2;
2816 for (nd = 1; nd <= i__1; ++nd) {
2818 for (ii = 1; ii <= i__2; ++ii) {
2821 for (kk = 1; kk <= i__3; ++kk) {
2822 bid1 = sotbv1[kk + (nd + ii * sotbv1_dim2) * sotbv1_dim1];
2823 bid2 = sotbv2[kk + (nd + ii * sotbv2_dim2) * sotbv2_dim1];
2824 bid3 = ditbv1[kk + (nd + ii * ditbv1_dim2) * ditbv1_dim1];
2825 bid4 = ditbv2[kk + (nd + ii * ditbv2_dim2) * ditbv2_dim1];
2827 for (jj = 1; jj <= i__4; ++jj) {
2828 jjp = (*nbpntv + 1) / 2 + jj;
2829 jjm = nvroo - jj + 1;
2830 sosotb[kk + (jj + nd * sosotb_dim2) * sosotb_dim1] =
2831 sosotb[kk + (jj + nd * sosotb_dim2) * sosotb_dim1]
2832 - bid1 * (fpntab[jjp + ((ii << 1) - 1) *
2833 fpntab_dim1] + fpntab[jjm + ((ii << 1) - 1) *
2834 fpntab_dim1]) - bid2 * (fpntab[jjp + (ii << 1) *
2835 fpntab_dim1] + fpntab[jjm + (ii << 1) *
2837 disotb[kk + (jj + nd * disotb_dim2) * disotb_dim1] =
2838 disotb[kk + (jj + nd * disotb_dim2) * disotb_dim1]
2839 - bid3 * (fpntab[jjp + ((ii << 1) - 1) *
2840 fpntab_dim1] + fpntab[jjm + ((ii << 1) - 1) *
2841 fpntab_dim1]) - bid4 * (fpntab[jjp + (ii << 1) *
2842 fpntab_dim1] + fpntab[jjm + (ii << 1) *
2844 soditb[kk + (jj + nd * soditb_dim2) * soditb_dim1] =
2845 soditb[kk + (jj + nd * soditb_dim2) * soditb_dim1]
2846 - bid1 * (fpntab[jjp + ((ii << 1) - 1) *
2847 fpntab_dim1] - fpntab[jjm + ((ii << 1) - 1) *
2848 fpntab_dim1]) - bid2 * (fpntab[jjp + (ii << 1) *
2849 fpntab_dim1] - fpntab[jjm + (ii << 1) *
2851 diditb[kk + (jj + nd * diditb_dim2) * diditb_dim1] =
2852 diditb[kk + (jj + nd * diditb_dim2) * diditb_dim1]
2853 - bid3 * (fpntab[jjp + ((ii << 1) - 1) *
2854 fpntab_dim1] - fpntab[jjm + ((ii << 1) - 1) *
2855 fpntab_dim1]) - bid4 * (fpntab[jjp + (ii << 1) *
2856 fpntab_dim1] - fpntab[jjm + (ii << 1) *
2865 /* ------------ Case when the discretization is done only on the roots */
2866 /* ---------- of Legendre polynom of uneven degree, 0 is root */
2869 if (*nbpntv % 2 == 1) {
2871 for (ii = 1; ii <= i__2; ++ii) {
2873 for (kk = 1; kk <= i__3; ++kk) {
2874 bid1 = sotbv1[kk + (nd + ii * sotbv1_dim2) * sotbv1_dim1]
2875 * fpntab[nvroo + 1 + ((ii << 1) - 1) *
2876 fpntab_dim1] + sotbv2[kk + (nd + ii * sotbv2_dim2)
2877 * sotbv2_dim1] * fpntab[nvroo + 1 + (ii << 1) *
2879 sosotb[kk + nd * sosotb_dim2 * sosotb_dim1] -= bid1;
2880 bid2 = ditbv1[kk + (nd + ii * ditbv1_dim2) * ditbv1_dim1]
2881 * fpntab[nvroo + 1 + ((ii << 1) - 1) *
2882 fpntab_dim1] + ditbv2[kk + (nd + ii * ditbv2_dim2)
2883 * ditbv2_dim1] * fpntab[nvroo + 1 + (ii << 1) *
2885 diditb[kk + nd * diditb_dim2 * diditb_dim1] -= bid2;
2892 if (*nbpntu % 2 == 1) {
2894 for (ii = 1; ii <= i__2; ++ii) {
2896 for (jj = 1; jj <= i__3; ++jj) {
2897 jjp = (*nbpntv + 1) / 2 + jj;
2898 jjm = nvroo - jj + 1;
2899 bid1 = sotbv1[(nd + ii * sotbv1_dim2) * sotbv1_dim1] * (
2900 fpntab[jjp + ((ii << 1) - 1) * fpntab_dim1] +
2901 fpntab[jjm + ((ii << 1) - 1) * fpntab_dim1]) +
2902 sotbv2[(nd + ii * sotbv2_dim2) * sotbv2_dim1] * (
2903 fpntab[jjp + (ii << 1) * fpntab_dim1] + fpntab[
2904 jjm + (ii << 1) * fpntab_dim1]);
2905 sosotb[(jj + nd * sosotb_dim2) * sosotb_dim1] -= bid1;
2906 bid2 = sotbv1[(nd + ii * sotbv1_dim2) * sotbv1_dim1] * (
2907 fpntab[jjp + ((ii << 1) - 1) * fpntab_dim1] -
2908 fpntab[jjm + ((ii << 1) - 1) * fpntab_dim1]) +
2909 sotbv2[(nd + ii * sotbv2_dim2) * sotbv2_dim1] * (
2910 fpntab[jjp + (ii << 1) * fpntab_dim1] - fpntab[
2911 jjm + (ii << 1) * fpntab_dim1]);
2912 diditb[jj + nd * diditb_dim2 * diditb_dim1] -= bid2;
2919 if (*nbpntu % 2 == 1 && *nbpntv % 2 == 1) {
2921 for (ii = 1; ii <= i__2; ++ii) {
2922 bid1 = sotbv1[(nd + ii * sotbv1_dim2) * sotbv1_dim1] * fpntab[
2923 nvroo + 1 + ((ii << 1) - 1) * fpntab_dim1] + sotbv2[(
2924 nd + ii * sotbv2_dim2) * sotbv2_dim1] * fpntab[nvroo
2925 + 1 + (ii << 1) * fpntab_dim1];
2926 sosotb[nd * sosotb_dim2 * sosotb_dim1] -= bid1;
2935 /* ------------------------------ The End -------------------------------
2940 AdvApp2Var_SysBase::mgsomsg_("MMA2CD2", 7L);
2945 //=======================================================================
2946 //function : mma2cd3_
2948 //=======================================================================
2949 int mma2cd3_(integer *ndimen,
2966 static integer c__1 = 1;
2968 /* System generated locals */
2969 integer sotbu1_dim1, sotbu1_dim2, sotbu1_offset, sotbu2_dim1, sotbu2_dim2,
2970 sotbu2_offset, ditbu1_dim1, ditbu1_dim2, ditbu1_offset,
2971 ditbu2_dim1, ditbu2_dim2, ditbu2_offset, fpntab_dim1,
2972 fpntab_offset, uhermt_dim1, uhermt_offset, sosotb_dim1,
2973 sosotb_dim2, sosotb_offset, diditb_dim1, diditb_dim2,
2974 diditb_offset, soditb_dim1, soditb_dim2, soditb_offset,
2975 disotb_dim1, disotb_dim2, disotb_offset, i__1, i__2, i__3, i__4;
2977 /* Local variables */
2978 static integer ncfhu, nuroo, nvroo, ii, nd, jj, kk, ibb, kkm, kkp;
2979 static doublereal bid1, bid2, bid3, bid4;
2981 /* **********************************************************************
2985 /* Discretisation on the parameters of polynoms of interpolation */
2986 /* of constraints on 2 borders iso-U of order IORDRU. */
2991 /* TOUS, AB_SPECIFI::CONTRAINTE&, DISCRETISATION, &POINT */
2993 /* INPUT ARGUMENTS : */
2994 /* ------------------ */
2995 /* NDIMEN: Dimension of the space. */
2996 /* NBPNTU: Nb of INTERNAL parameters of discretisation by U. */
2997 /* This is also the nb of root of Legendre polynom where discretization is done. */
2998 /* UROOTL: Table of parameters of discretisation ON (-1,1) by U.
3000 /* NBPNTV: Nb of INTERNAL parameters of discretisation by V. */
3001 /* This is also the nb of root of Legendre polynom where discretization is done. */
3002 /* IORDRV: Order of constraint imposed at the extremities of iso-V */
3003 /* = 0, calculate the extremities of iso-V */
3004 /* = 1, calculate, additionally, the 1st derivative in the direction of iso-V */
3005 /* = 2, calculate, additionally, the 2nd derivative in the direction of iso-V */
3006 /* SOTBU1: Table of NBPNTU/2 sums of 2 index points */
3007 /* NBPNTU-II+1 and II, for II = 1, NBPNTU/2 on iso-V0. */
3008 /* SOTBU2: Table of NBPNTV/2 sums of 2 index points */
3009 /* NBPNTU-II+1 and II, for II = 1, NBPNTU/2 on iso-V1. */
3010 /* DITBU1: Table of NBPNTU/2 differences of 2 index points */
3011 /* NBPNTU-II+1 and II, for II = 1, NBPNTU/2 on iso-V0. */
3012 /* DITBU2: Table of NBPNTU/2 differences of 2 index points */
3013 /* NBPNTU-II+1 and II, for II = 1, NBPNTU/2 on iso-V1. */
3014 /* SOSOTB: Preinitialized table (input/output argument). */
3015 /* DISOTB: Preinitialized table (input/output argument). */
3016 /* SODITB: Preinitialized table (input/output argument). */
3017 /* DIDITB: Preinitialized table (input/output argument) */
3019 /* OUTPUT ARGUMENTS : */
3020 /* ------------------- */
3021 /* FPNTAB: Auxiliary table. */
3022 /* UHERMT: Table of 2*(IORDRU+1) coeff. of 2*(IORDRU+1) polynoms of Hermite. */
3023 /* SOSOTB: Table where the terms of constraints are added */
3024 /* C(ui,vj) + C(ui,-vj) + C(-ui,vj) + C(-ui,-vj) */
3025 /* with ui and vj positive roots of the Legendre polynom */
3026 /* of degree NBPNTU and NBPNTV respectively. */
3027 /* DISOTB: Table where the terms of constraints are added */
3028 /* C(ui,vj) + C(ui,-vj) - C(-ui,vj) - C(-ui,-vj) */
3029 /* with ui and vj positive roots of the polynom of Legendre */
3030 /* of degree NBPNTU and NBPNTV respectively. */
3031 /* SODITB: Table where the terms of constraints are added */
3032 /* C(ui,vj) - C(ui,-vj) + C(-ui,vj) - C(-ui,-vj) */
3033 /* with ui and vj positive roots of the polynom of Legendre */
3034 /* of degree NBPNTU and NBPNTV respectively. */
3035 /* DIDITB: Table where the terms of constraints are added */
3036 /* C(ui,vj) - C(ui,-vj) - C(-ui,vj) + C(-ui,-vj) */
3037 /* with ui and vj positive roots of the polynom of Legendre */
3038 /* of degree NBPNTU and NBPNTV respectively. */
3040 /* COMMONS USED : */
3041 /* ---------------- */
3043 /* REFERENCES CALLED : */
3044 /* ----------------------- */
3046 /* DESCRIPTION/NOTES/LIMITATIONS : */
3047 /* ----------------------------------- */
3049 /* $ HISTORIQUE DES MODIFICATIONS : */
3050 /* -------------------------------- */
3051 /* 08-08-1991: RBD; Creation. */
3053 /* **********************************************************************
3056 /* Name of the routine */
3059 /* Parameter adjustments */
3061 diditb_dim1 = *nbpntu / 2 + 1;
3062 diditb_dim2 = *nbpntv / 2 + 1;
3063 diditb_offset = diditb_dim1 * diditb_dim2;
3064 diditb -= diditb_offset;
3065 disotb_dim1 = *nbpntu / 2;
3066 disotb_dim2 = *nbpntv / 2;
3067 disotb_offset = disotb_dim1 * (disotb_dim2 + 1) + 1;
3068 disotb -= disotb_offset;
3069 soditb_dim1 = *nbpntu / 2;
3070 soditb_dim2 = *nbpntv / 2;
3071 soditb_offset = soditb_dim1 * (soditb_dim2 + 1) + 1;
3072 soditb -= soditb_offset;
3073 sosotb_dim1 = *nbpntu / 2 + 1;
3074 sosotb_dim2 = *nbpntv / 2 + 1;
3075 sosotb_offset = sosotb_dim1 * sosotb_dim2;
3076 sosotb -= sosotb_offset;
3077 uhermt_dim1 = (*iordru << 1) + 2;
3078 uhermt_offset = uhermt_dim1;
3079 uhermt -= uhermt_offset;
3080 fpntab_dim1 = *nbpntu;
3081 fpntab_offset = fpntab_dim1 + 1;
3082 fpntab -= fpntab_offset;
3083 ditbu2_dim1 = *nbpntv / 2 + 1;
3084 ditbu2_dim2 = *ndimen;
3085 ditbu2_offset = ditbu2_dim1 * (ditbu2_dim2 + 1);
3086 ditbu2 -= ditbu2_offset;
3087 ditbu1_dim1 = *nbpntv / 2 + 1;
3088 ditbu1_dim2 = *ndimen;
3089 ditbu1_offset = ditbu1_dim1 * (ditbu1_dim2 + 1);
3090 ditbu1 -= ditbu1_offset;
3091 sotbu2_dim1 = *nbpntv / 2 + 1;
3092 sotbu2_dim2 = *ndimen;
3093 sotbu2_offset = sotbu2_dim1 * (sotbu2_dim2 + 1);
3094 sotbu2 -= sotbu2_offset;
3095 sotbu1_dim1 = *nbpntv / 2 + 1;
3096 sotbu1_dim2 = *ndimen;
3097 sotbu1_offset = sotbu1_dim1 * (sotbu1_dim2 + 1);
3098 sotbu1 -= sotbu1_offset;
3101 ibb = AdvApp2Var_SysBase::mnfndeb_();
3103 AdvApp2Var_SysBase::mgenmsg_("MMA2CD3", 7L);
3106 /* ------------------- Discretization of polynoms of Hermit -----------
3109 ncfhu = (*iordru + 1) << 1;
3111 for (ii = 1; ii <= i__1; ++ii) {
3113 for (kk = 1; kk <= i__2; ++kk) {
3114 AdvApp2Var_MathBase::mmmpocur_(&ncfhu,
3117 &uhermt[ii * uhermt_dim1],
3119 &fpntab[kk + ii * fpntab_dim1]);
3125 /* ---- The discretizations of polynoms of constraints are subtracted ----
3128 nvroo = *nbpntv / 2;
3129 nuroo = *nbpntu / 2;
3132 for (nd = 1; nd <= i__1; ++nd) {
3134 for (ii = 1; ii <= i__2; ++ii) {
3137 for (jj = 1; jj <= i__3; ++jj) {
3138 bid1 = sotbu1[jj + (nd + ii * sotbu1_dim2) * sotbu1_dim1];
3139 bid2 = sotbu2[jj + (nd + ii * sotbu2_dim2) * sotbu2_dim1];
3140 bid3 = ditbu1[jj + (nd + ii * ditbu1_dim2) * ditbu1_dim1];
3141 bid4 = ditbu2[jj + (nd + ii * ditbu2_dim2) * ditbu2_dim1];
3143 for (kk = 1; kk <= i__4; ++kk) {
3144 kkp = (*nbpntu + 1) / 2 + kk;
3145 kkm = nuroo - kk + 1;
3146 sosotb[kk + (jj + nd * sosotb_dim2) * sosotb_dim1] =
3147 sosotb[kk + (jj + nd * sosotb_dim2) * sosotb_dim1]
3148 - bid1 * (fpntab[kkp + ((ii << 1) - 1) *
3149 fpntab_dim1] + fpntab[kkm + ((ii << 1) - 1) *
3150 fpntab_dim1]) - bid2 * (fpntab[kkp + (ii << 1) *
3151 fpntab_dim1] + fpntab[kkm + (ii << 1) *
3153 disotb[kk + (jj + nd * disotb_dim2) * disotb_dim1] =
3154 disotb[kk + (jj + nd * disotb_dim2) * disotb_dim1]
3155 - bid1 * (fpntab[kkp + ((ii << 1) - 1) *
3156 fpntab_dim1] - fpntab[kkm + ((ii << 1) - 1) *
3157 fpntab_dim1]) - bid2 * (fpntab[kkp + (ii << 1) *
3158 fpntab_dim1] - fpntab[kkm + (ii << 1) *
3160 soditb[kk + (jj + nd * soditb_dim2) * soditb_dim1] =
3161 soditb[kk + (jj + nd * soditb_dim2) * soditb_dim1]
3162 - bid3 * (fpntab[kkp + ((ii << 1) - 1) *
3163 fpntab_dim1] + fpntab[kkm + ((ii << 1) - 1) *
3164 fpntab_dim1]) - bid4 * (fpntab[kkp + (ii << 1) *
3165 fpntab_dim1] + fpntab[kkm + (ii << 1) *
3167 diditb[kk + (jj + nd * diditb_dim2) * diditb_dim1] =
3168 diditb[kk + (jj + nd * diditb_dim2) * diditb_dim1]
3169 - bid3 * (fpntab[kkp + ((ii << 1) - 1) *
3170 fpntab_dim1] - fpntab[kkm + ((ii << 1) - 1) *
3171 fpntab_dim1]) - bid4 * (fpntab[kkp + (ii << 1) *
3172 fpntab_dim1] - fpntab[kkm + (ii << 1) *
3181 /* ------------ Case when the discretization is done only on the roots */
3182 /* ---------- of Legendre polynom of uneven degree, 0 is root */
3186 if (*nbpntu % 2 == 1) {
3188 for (ii = 1; ii <= i__2; ++ii) {
3190 for (jj = 1; jj <= i__3; ++jj) {
3191 bid1 = sotbu1[jj + (nd + ii * sotbu1_dim2) * sotbu1_dim1]
3192 * fpntab[nuroo + 1 + ((ii << 1) - 1) *
3193 fpntab_dim1] + sotbu2[jj + (nd + ii * sotbu2_dim2)
3194 * sotbu2_dim1] * fpntab[nuroo + 1 + (ii << 1) *
3196 sosotb[(jj + nd * sosotb_dim2) * sosotb_dim1] -= bid1;
3197 bid2 = ditbu1[jj + (nd + ii * ditbu1_dim2) * ditbu1_dim1]
3198 * fpntab[nuroo + 1 + ((ii << 1) - 1) *
3199 fpntab_dim1] + ditbu2[jj + (nd + ii * ditbu2_dim2)
3200 * ditbu2_dim1] * fpntab[nuroo + 1 + (ii << 1) *
3202 diditb[(jj + nd * diditb_dim2) * diditb_dim1] -= bid2;
3209 if (*nbpntv % 2 == 1) {
3211 for (ii = 1; ii <= i__2; ++ii) {
3213 for (kk = 1; kk <= i__3; ++kk) {
3214 kkp = (*nbpntu + 1) / 2 + kk;
3215 kkm = nuroo - kk + 1;
3216 bid1 = sotbu1[(nd + ii * sotbu1_dim2) * sotbu1_dim1] * (
3217 fpntab[kkp + ((ii << 1) - 1) * fpntab_dim1] +
3218 fpntab[kkm + ((ii << 1) - 1) * fpntab_dim1]) +
3219 sotbu2[(nd + ii * sotbu2_dim2) * sotbu2_dim1] * (
3220 fpntab[kkp + (ii << 1) * fpntab_dim1] + fpntab[
3221 kkm + (ii << 1) * fpntab_dim1]);
3222 sosotb[kk + nd * sosotb_dim2 * sosotb_dim1] -= bid1;
3223 bid2 = sotbu1[(nd + ii * sotbu1_dim2) * sotbu1_dim1] * (
3224 fpntab[kkp + ((ii << 1) - 1) * fpntab_dim1] -
3225 fpntab[kkm + ((ii << 1) - 1) * fpntab_dim1]) +
3226 sotbu2[(nd + ii * sotbu2_dim2) * sotbu2_dim1] * (
3227 fpntab[kkp + (ii << 1) * fpntab_dim1] - fpntab[
3228 kkm + (ii << 1) * fpntab_dim1]);
3229 diditb[kk + nd * diditb_dim2 * diditb_dim1] -= bid2;
3236 if (*nbpntu % 2 == 1 && *nbpntv % 2 == 1) {
3238 for (ii = 1; ii <= i__2; ++ii) {
3239 bid1 = sotbu1[(nd + ii * sotbu1_dim2) * sotbu1_dim1] * fpntab[
3240 nuroo + 1 + ((ii << 1) - 1) * fpntab_dim1] + sotbu2[(
3241 nd + ii * sotbu2_dim2) * sotbu2_dim1] * fpntab[nuroo
3242 + 1 + (ii << 1) * fpntab_dim1];
3243 sosotb[nd * sosotb_dim2 * sosotb_dim1] -= bid1;
3252 /* ------------------------------ The End -------------------------------
3257 AdvApp2Var_SysBase::mgsomsg_("MMA2CD3", 7L);
3262 //=======================================================================
3263 //function : mma2cdi_
3265 //=======================================================================
3266 int AdvApp2Var_ApproxF2var::mma2cdi_( integer *ndimen,
3292 static integer c__8 = 8;
3294 /* System generated locals */
3295 integer contr1_dim1, contr1_dim2, contr1_offset, contr2_dim1, contr2_dim2,
3296 contr2_offset, contr3_dim1, contr3_dim2, contr3_offset,
3297 contr4_dim1, contr4_dim2, contr4_offset, sosotb_dim1, sosotb_dim2,
3298 sosotb_offset, diditb_dim1, diditb_dim2, diditb_offset,
3299 soditb_dim1, soditb_dim2, soditb_offset, disotb_dim1, disotb_dim2,
3302 /* Local variables */
3303 static integer ilong;
3304 static long int iofwr;
3305 static doublereal wrkar[1];
3306 static integer iszwr;
3307 static integer ibb, ier;
3308 static integer isz1, isz2, isz3, isz4;
3309 static long int ipt1, ipt2, ipt3, ipt4;
3314 /* **********************************************************************
3319 /* Discretisation on the parameters of polynomes of interpolation */
3320 /* of constraints of order IORDRE. */
3324 /* TOUS, AB_SPECIFI::CONTRAINTE&, DISCRETISATION, &POINT */
3326 //* INPUT ARGUMENTS : */
3327 /* ------------------ */
3328 /* NDIMEN: Dimension of the space. */
3329 /* NBPNTU: Nb of INTERNAL parameters of discretisation by U. */
3330 /* This is also the nb of root of Legendre polynom where discretization is done. */
3331 /* UROOTL: Table of parameters of discretisation ON (-1,1) by U.
3333 /* NBPNTV: Nb of INTERNAL parameters of discretisation by V. */
3334 /* This is also the nb of root of Legendre polynom where discretization is done. */
3335 /* VROOTL: Table of parameters of discretisation ON (-1,1) by V.
3337 /* IORDRV: Order of constraint imposed at the extremities of iso-U */
3338 /* = 0, calculate the extremities of iso-U */
3339 /* = 1, calculate, additionally, the 1st derivative in the direction of iso-U */
3340 /* = 2, calculate, additionally, the 2nd derivative in the direction of iso-U */
3341 /* IORDRU: Order of constraint imposed at the extremities of iso-V */
3342 /* = 0, calculate the extremities of iso-V */
3343 /* = 1, calculate, additionally, the 1st derivative in the direction of iso-V */
3344 /* = 2, calculate, additionally, the 2nd derivative in the direction of iso-V */
3345 /* CONTR1: Contains, if IORDRU and IORDRV>=0, the values at the */
3346 /* extremities of F(U0,V0) and its derivatives. */
3347 /* CONTR2: Contains, if IORDRU and IORDRV>=0, the values at the */
3348 /* extremities of F(U1,V0) and its derivatives. */
3349 /* CONTR3: Contains, if IORDRU and IORDRV>=0, the values at the */
3350 /* extremities of F(U0,V1) and its derivatives. */
3351 /* CONTR4: Contains, if IORDRU and IORDRV>=0, the values at the */
3352 /* extremities of F(U1,V1) and its derivatives. */
3353 /* SOTBU1: Table of NBPNTU/2 sums of 2 index points */
3354 /* NBPNTU-II+1 and II, for II = 1, NBPNTU/2 on iso-V0. */
3355 /* SOTBU2: Table of NBPNTV/2 sums of 2 index points */
3356 /* NBPNTU-II+1 and II, for II = 1, NBPNTU/2 on iso-V1. */
3357 /* DITBU1: Table of NBPNTU/2 differences of 2 index points */
3358 /* NBPNTU-II+1 and II, for II = 1, NBPNTU/2 on iso-V0. */
3359 /* DITBU2: Table of NBPNTU/2 differences of 2 index points */
3360 /* NBPNTU-II+1 and II, for II = 1, NBPNTU/2 on iso-V1. */
3361 /* SOTBV1: Table of NBPNTV/2 sums of 2 index points */
3362 /* NBPNTV-II+1 and II, for II = 1, NBPNTV/2 on iso-V0. */
3363 /* SOTBV2: Table of NBPNTV/2 sums of 2 index points */
3364 /* NBPNTV-II+1 and II, for II = 1, NBPNTV/2 on iso-V1. */
3365 /* DITBV1: Table of NBPNTV/2 differences of 2 index points */
3366 /* NBPNTV-II+1 and II, for II = 1, NBPNTV/2 on iso-V0. */
3367 /* DITBV2: Table of NBPNTV/2 differences of 2 index points */
3368 /* NBPNTV-II+1 and II, for II = 1, NBPNTV/2 on iso-V1. */
3369 /* SOSOTB: Preinitialized table (input/output argument). */
3370 /* DISOTB: Preinitialized table (input/output argument). */
3371 /* SODITB: Preinitialized table (input/output argument). */
3372 /* DIDITB: Preinitialized table (input/output argument) */
3374 /* ARGUMENTS DE SORTIE : */
3375 /* ------------------- */
3376 /* SOSOTB: Table where the terms of constraints are added */
3377 /* C(ui,vj) + C(ui,-vj) + C(-ui,vj) + C(-ui,-vj) */
3378 /* with ui and vj positive roots of the Legendre polynom */
3379 /* of degree NBPNTU and NBPNTV respectively. */
3380 /* DISOTB: Table where the terms of constraints are added */
3381 /* C(ui,vj) + C(ui,-vj) - C(-ui,vj) - C(-ui,-vj) */
3382 /* with ui and vj positive roots of the polynom of Legendre */
3383 /* of degree NBPNTU and NBPNTV respectively. */
3384 /* SODITB: Table where the terms of constraints are added */
3385 /* C(ui,vj) - C(ui,-vj) + C(-ui,vj) - C(-ui,-vj) */
3386 /* with ui and vj positive roots of the polynom of Legendre */
3387 /* of degree NBPNTU and NBPNTV respectively. */
3388 /* DIDITB: Table where the terms of constraints are added */
3389 /* C(ui,vj) - C(ui,-vj) - C(-ui,vj) + C(-ui,-vj) */
3390 /* with ui and vj positive roots of the polynom of Legendre */
3391 /* of degree NBPNTU and NBPNTV respectively. */
3392 /* IERCOD: = 0, OK, */
3393 /* = 1, Value or IORDRV or IORDRU is out of allowed values. */
3394 /* =13, Pb of dynamic allocation. */
3396 /* COMMONS USED : */
3397 /* ---------------- */
3399 /* REFERENCES CALLED : */
3400 /* -------------------- */
3402 /* DESCRIPTION/NOTES/LIMITATIONS : */
3403 /* ------------------------------- */
3406 /* **********************************************************************
3409 /* The name of the routine */
3412 /* Parameter adjustments */
3414 diditb_dim1 = *nbpntu / 2 + 1;
3415 diditb_dim2 = *nbpntv / 2 + 1;
3416 diditb_offset = diditb_dim1 * diditb_dim2;
3417 diditb -= diditb_offset;
3418 disotb_dim1 = *nbpntu / 2;
3419 disotb_dim2 = *nbpntv / 2;
3420 disotb_offset = disotb_dim1 * (disotb_dim2 + 1) + 1;
3421 disotb -= disotb_offset;
3422 soditb_dim1 = *nbpntu / 2;
3423 soditb_dim2 = *nbpntv / 2;
3424 soditb_offset = soditb_dim1 * (soditb_dim2 + 1) + 1;
3425 soditb -= soditb_offset;
3426 sosotb_dim1 = *nbpntu / 2 + 1;
3427 sosotb_dim2 = *nbpntv / 2 + 1;
3428 sosotb_offset = sosotb_dim1 * sosotb_dim2;
3429 sosotb -= sosotb_offset;
3431 contr4_dim1 = *ndimen;
3432 contr4_dim2 = *iordru + 2;
3433 contr4_offset = contr4_dim1 * (contr4_dim2 + 1) + 1;
3434 contr4 -= contr4_offset;
3435 contr3_dim1 = *ndimen;
3436 contr3_dim2 = *iordru + 2;
3437 contr3_offset = contr3_dim1 * (contr3_dim2 + 1) + 1;
3438 contr3 -= contr3_offset;
3439 contr2_dim1 = *ndimen;
3440 contr2_dim2 = *iordru + 2;
3441 contr2_offset = contr2_dim1 * (contr2_dim2 + 1) + 1;
3442 contr2 -= contr2_offset;
3443 contr1_dim1 = *ndimen;
3444 contr1_dim2 = *iordru + 2;
3445 contr1_offset = contr1_dim1 * (contr1_dim2 + 1) + 1;
3446 contr1 -= contr1_offset;
3457 ibb = AdvApp2Var_SysBase::mnfndeb_();
3459 AdvApp2Var_SysBase::mgenmsg_("MMA2CDI", 7L);
3463 if (*iordru < -1 || *iordru > 2) {
3466 if (*iordrv < -1 || *iordrv > 2) {
3470 /* ------------------------- Set to zero --------------------------------
3473 ilong = (*nbpntu / 2 + 1) * (*nbpntv / 2 + 1) * *ndimen;
3474 AdvApp2Var_SysBase::mvriraz_(&ilong, (char *)&sosotb[sosotb_offset]);
3475 AdvApp2Var_SysBase::mvriraz_(&ilong, (char *)&diditb[diditb_offset]);
3476 ilong = *nbpntu / 2 * (*nbpntv / 2) * *ndimen;
3477 AdvApp2Var_SysBase::mvriraz_(&ilong, (char *)&soditb[soditb_offset]);
3478 AdvApp2Var_SysBase::mvriraz_(&ilong, (char *)&disotb[disotb_offset]);
3479 if (*iordru == -1 && *iordrv == -1) {
3485 isz1 = ((*iordru + 1) << 2) * (*iordru + 1);
3486 isz2 = ((*iordrv + 1) << 2) * (*iordrv + 1);
3487 isz3 = ((*iordru + 1) << 1) * *nbpntu;
3488 isz4 = ((*iordrv + 1) << 1) * *nbpntv;
3489 iszwr = isz1 + isz2 + isz3 + isz4;
3490 AdvApp2Var_SysBase::mcrrqst_(&c__8, &iszwr, wrkar, &iofwr, &ier);
3499 if (*iordru >= 0 && *iordru <= 2) {
3501 /* --- Return 2*(IORDRU+1) coeff of 2*(IORDRU+1) polynoms of Hermite
3504 AdvApp2Var_ApproxF2var::mma1her_(iordru, &wrkar[ipt1], iercod);
3509 /* ---- Subract discretizations of polynoms of constraints
3512 mma2cd3_(ndimen, nbpntu, &urootl[1], nbpntv, iordru, &sotbu1[1], &
3513 sotbu2[1], &ditbu1[1], &ditbu2[1], &wrkar[ipt3], &wrkar[ipt1],
3514 &sosotb[sosotb_offset], &soditb[soditb_offset], &disotb[
3515 disotb_offset], &diditb[diditb_offset]);
3518 if (*iordrv >= 0 && *iordrv <= 2) {
3520 /* --- Return 2*(IORDRV+1) coeff of 2*(IORDRV+1) polynoms of Hermite
3523 AdvApp2Var_ApproxF2var::mma1her_(iordrv, &wrkar[ipt2], iercod);
3528 /* ---- Subtract discretisations of polynoms of constraint
3531 mma2cd2_(ndimen, nbpntu, nbpntv, &vrootl[1], iordrv, &sotbv1[1], &
3532 sotbv2[1], &ditbv1[1], &ditbv2[1], &wrkar[ipt4], &wrkar[ipt2],
3533 &sosotb[sosotb_offset], &soditb[soditb_offset], &disotb[
3534 disotb_offset], &diditb[diditb_offset]);
3537 /* --------------- Subtract constraints of corners ----------------
3540 if (*iordru >= 0 && *iordrv >= 0) {
3541 mma2cd1_(ndimen, nbpntu, &urootl[1], nbpntv, &vrootl[1], iordru,
3542 iordrv, &contr1[contr1_offset], &contr2[contr2_offset], &
3543 contr3[contr3_offset], &contr4[contr4_offset], &wrkar[ipt3], &
3544 wrkar[ipt4], &wrkar[ipt1], &wrkar[ipt2], &sosotb[
3545 sosotb_offset], &soditb[soditb_offset], &disotb[disotb_offset]
3546 , &diditb[diditb_offset]);
3550 /* ------------------------------ The End -------------------------------
3552 /* --> IORDRE is not within the autorised diapason. */
3556 /* --> PB of dynamic allocation. */
3563 AdvApp2Var_SysBase::mcrdelt_(&c__8, &iszwr, wrkar, &iofwr, &ier);
3568 AdvApp2Var_SysBase::maermsg_("MMA2CDI", iercod, 7L);
3570 AdvApp2Var_SysBase::mgsomsg_("MMA2CDI", 7L);
3575 //=======================================================================
3576 //function : mma2ce1_
3578 //=======================================================================
3579 int AdvApp2Var_ApproxF2var::mma2ce1_(integer *numdec,
3607 static integer c__8 = 8;
3609 /* System generated locals */
3610 integer sosotb_dim1, sosotb_dim2, sosotb_offset, disotb_dim1, disotb_dim2,
3611 disotb_offset, soditb_dim1, soditb_dim2, soditb_offset,
3612 diditb_dim1, diditb_dim2, diditb_offset, patjac_dim1, patjac_dim2,
3615 /* Local variables */
3616 static logical ldbg;
3617 static long int iofwr;
3618 static doublereal wrkar[1];
3619 static integer iszwr;
3621 static integer isz1, isz2, isz3, isz4, isz5, isz6, isz7;
3622 static long int ipt1, ipt2, ipt3, ipt4, ipt5, ipt6, ipt7;
3626 /* **********************************************************************
3631 /* Calculation of coefficients of polynomial approximation of degree */
3632 /* (NDJACU,NDJACV) of a function F(u,v), starting from its */
3633 /* discretization on roots of Legendre polynom of degree */
3634 /* NBPNTU by U and NBPNTV by V. */
3638 /* TOUS,AB_SPECIFI::FONCTION&,APPROXIMATION,&POLYNOME,&ERREUR */
3640 /* INPUT ARGUMENTS : */
3641 /* ------------------ */
3642 /* NUMDEC: Indicates if it is POSSIBLE to cut function F(u,v). */
3643 /* = 5, It is POSSIBLE to cut by U or by V or in both directions simultaneously. */
3644 /* = 4, It is POSSIBLE to cut by U or by V BUT NOT in both */
3645 /* directions simultaneously (cutting by V is preferable). */
3646 /* = 3, It is POSSIBLE to cut by U or by V BUT NOT in both */
3647 /* directions simultaneously (cutting by U is preferable). */
3648 /* = 2, It is POSSIBLE to cut only by V (i.e. insert parameter */
3649 /* of cutting Vj). */
3650 /* = 1, It is POSSIBLE to cut only by U (i.e. insert parameter */
3651 /* of cutting Ui). */
3652 /* = 0, It is not POSSIBLE to cut anything */
3653 /* NDIMEN: Dimension of the space. */
3654 /* NBSESP: Nb of independent sub-spaces on which the errors are calculated. */
3655 /* NDIMSE: Table of dimensions of each of sub-spaces. */
3656 /* NDMINU: Minimum degree by U to be preserved for the approximation. */
3657 /* NDMINV: Minimum degree by V to be preserved for the approximation. */
3658 /* NDGULI: Limit of nb of coefficients by U of the solution. */
3659 /* NDGVLI: Limit of nb of coefficients by V of the solution. */
3660 /* NDJACU: Max degree of the polynom of approximation by U. */
3661 /* The representation in the orthogonal base starts from degree */
3662 /* 0 to degree NDJACU-2*(IORDRU+1). The polynomial base is the base of */
3663 /* Jacobi of order -1 (Legendre), 0, 1 or 2. */
3664 /* It is required that 2*IORDRU+1 <= NDMINU <= NDGULI < NDJACU */
3665 /* NDJACV: Max degree of the polynom of approximation by V. */
3666 /* The representation in the orthogonal base starts from degree */
3667 /* 0 to degree NDJACV-2*(IORDRV+1). The polynomial base is */
3668 /* the base of Jacobi of order -1 (Legendre), 0, 1 or 2 */
3669 /* It is required that 2*IORDRV+1 <= NDMINV <= NDGVLI < NDJACV */
3670 /* IORDRU: Order of the Jacobi base (-1,0,1 or 2) by U. Corresponds */
3671 /* to the step of constraints C0, C1 or C2. */
3672 /* IORDRV: Order of the Jacobi base (-1,0,1 or 2) by U. Corresponds */
3673 /* to the step of constraints C0, C1 or C2. */
3674 /* NBPNTU: Degree of Legendre polynom on the roots which of are */
3675 /* calculated the coefficients of integration by u */
3676 /* by Gauss method. It is required that NBPNTU = 30, 40, */
3677 /* 50 or 61 and NDJACU-2*(IORDRU+1) < NBPNTU. */
3678 /* NBPNTV: Degree of Legendre polynom on the roots which of are */
3679 /* calculated the coefficients of integration by u */
3680 /* by Gauss method. It is required that NBPNTV = 30, 40, */
3681 /* 50 or 61 and NDJACV-2*(IORDRV+1) < NBPNTV. */
3682 /* EPSAPR: Table of NBSESP tolerances imposed on each sub-spaces. */
3683 /* SOSOTB: Table of F(ui,vj) + F(ui,-vj) + F(-ui,vj) + F(-ui,-vj) */
3684 /* with ui and vj - positive roots of the Legendre polynom */
3685 /* of degree NBPNTU and NBPNTV respectively. Additionally, */
3686 /* table SOSOTB(0,j) contains F(0,vj) + F(0,-vj), */
3687 /* table SOSOTB(i,0) contains F(ui,0) + F(-ui,0) and */
3688 /* SOSOTB(0,0) contains F(0,0). */
3689 /* DISOTB: Table of F(ui,vj) + F(ui,-vj) - F(-ui,vj) - F(-ui,-vj) */
3690 /* with ui and vj positive roots of Legendre polynom */
3691 /* of degree NBPNTU and NBPNTV respectively. */
3692 /* SODITB: Table of F(ui,vj) - F(ui,-vj) + F(-ui,vj) - F(-ui,-vj) */
3693 /* with ui and vj positive roots of Legendre polynom */
3694 /* of degree NBPNTU and NBPNTV respectively. */
3695 /* DIDITB: Table of F(ui,vj) - F(ui,-vj) - F(-ui,vj) + F(-ui,-vj) */
3696 /* with ui and vj positive roots of Legendre polynom */
3697 /* of degree NBPNTU and NBPNTV respectively. Additionally, */
3698 /* table DIDITB(0,j) contains F(0,vj) - F(0,-vj), */
3699 /* and table DIDITB(i,0) contains F(ui,0) - F(-ui,0). */
3701 /* OUTPUT ARGUMENTS */
3702 /* --------------- */
3703 /* PATJAC: Table of coefficients of polynom P(u,v) of approximation */
3704 /* of F(u,v) with eventually taking into account of */
3705 /* constraints. P(u,v) is of degree (NDJACU,NDJACV). */
3706 /* This table contains other coeff if ITYDEC = 0. */
3707 /* ERRMAX: For 1<=i<=NBSESP, ERRMAX(i) contains max errors */
3708 /* on each of sub-spaces SI ITYDEC = 0. */
3709 /* ERRMOY: Contains average errors for each of NBSESP sub-spaces SI ITYDEC = 0. */
3710 /* NDEGPU: Degree by U for square PATJAC. Valable if ITYDEC=0. */
3711 /* NDEGPV: Degree by V for square PATJAC. Valable if ITYDEC=0. */
3712 /* ITYDEC: Shows if it is NECESSARY to cut again function F(u,v). */
3713 /* = 0, it is not NECESSARY to cut anything, PATJAC is OK. */
3714 /* = 1, it is NECESSARY to cut only by U (i.e. insert parameter of cutting Ui). */
3715 /* = 2, it is NECESSARY to cut only by V (i.e. insert parameter of cutting Vj). */
3716 /* = 3, it is NECESSARY to cut both by U AND by V. */
3717 /* IERCOD: Error code. */
3718 /* = 0, Everything is OK. */
3719 /* = -1, There is the best possible solution, but the */
3720 /* user tolerance is not satisfactory (3*only) */
3721 /* = 1, Incoherent entries. */
3723 /* COMMONS USED : */
3724 /* ---------------- */
3726 /* REFERENCES CALLED : */
3727 /* --------------------- */
3729 /* DESCRIPTION/NOTES/LIMITATIONS : */
3730 /* ------------------------------- */
3733 /* **********************************************************************
3735 /* Name of the routine */
3738 /* --------------------------- Initialisations --------------------------
3741 /* Parameter adjustments */
3746 patjac_dim1 = *ndjacu + 1;
3747 patjac_dim2 = *ndjacv + 1;
3748 patjac_offset = patjac_dim1 * patjac_dim2;
3749 patjac -= patjac_offset;
3750 diditb_dim1 = *nbpntu / 2 + 1;
3751 diditb_dim2 = *nbpntv / 2 + 1;
3752 diditb_offset = diditb_dim1 * diditb_dim2;
3753 diditb -= diditb_offset;
3754 soditb_dim1 = *nbpntu / 2;
3755 soditb_dim2 = *nbpntv / 2;
3756 soditb_offset = soditb_dim1 * (soditb_dim2 + 1) + 1;
3757 soditb -= soditb_offset;
3758 disotb_dim1 = *nbpntu / 2;
3759 disotb_dim2 = *nbpntv / 2;
3760 disotb_offset = disotb_dim1 * (disotb_dim2 + 1) + 1;
3761 disotb -= disotb_offset;
3762 sosotb_dim1 = *nbpntu / 2 + 1;
3763 sosotb_dim2 = *nbpntv / 2 + 1;
3764 sosotb_offset = sosotb_dim1 * sosotb_dim2;
3765 sosotb -= sosotb_offset;
3768 ldbg = AdvApp2Var_SysBase::mnfndeb_() >= 3;
3770 AdvApp2Var_SysBase::mgenmsg_("MMA2CE1", 7L);
3775 isz1 = (*nbpntu / 2 + 1) * (*ndjacu - ((*iordru + 1) << 1) + 1);
3776 isz2 = (*nbpntv / 2 + 1) * (*ndjacv - ((*iordrv + 1) << 1) + 1);
3777 isz3 = (*nbpntv / 2 + 1) * (*ndjacu - ((*iordru + 1) << 1) + 1) * *ndimen;
3778 isz4 = *nbpntv / 2 * (*ndjacu - ((*iordru + 1) << 1) + 1) * *ndimen;
3779 isz5 = *ndjacu + 1 - ((*iordru + 1) << 1);
3780 isz6 = *ndjacv + 1 - ((*iordrv + 1) << 1);
3781 isz7 = *ndimen << 2;
3782 iszwr = isz1 + isz2 + isz3 + isz4 + isz5 + isz6 + isz7;
3783 AdvApp2Var_SysBase::mcrrqst_(&c__8, &iszwr, wrkar, &iofwr, &ier);
3795 /* ----------------- Return Gauss coefficients of integration ----------------
3798 AdvApp2Var_ApproxF2var::mmapptt_(ndjacu, nbpntu, iordru, &wrkar[ipt1], iercod);
3802 AdvApp2Var_ApproxF2var::mmapptt_(ndjacv, nbpntv, iordrv, &wrkar[ipt2], iercod);
3807 /* ------------------- Return max polynoms of Jacobi ------------
3810 AdvApp2Var_ApproxF2var::mma2jmx_(ndjacu, iordru, &wrkar[ipt5]);
3811 AdvApp2Var_ApproxF2var::mma2jmx_(ndjacv, iordrv, &wrkar[ipt6]);
3813 /* ------ Calculate the coefficients and their contribution to the error ----
3816 mma2ce2_(numdec, ndimen, nbsesp, &ndimse[1], ndminu, ndminv, ndguli,
3817 ndgvli, ndjacu, ndjacv, iordru, iordrv, nbpntu, nbpntv, &epsapr[1]
3818 , &sosotb[sosotb_offset], &disotb[disotb_offset], &soditb[
3819 soditb_offset], &diditb[diditb_offset], &wrkar[ipt1], &wrkar[ipt2]
3820 , &wrkar[ipt5], &wrkar[ipt6], &wrkar[ipt7], &wrkar[ipt3], &wrkar[
3821 ipt4], &patjac[patjac_offset], &errmax[1], &errmoy[1], ndegpu,
3822 ndegpv, itydec, iercod);
3828 /* ------------------------------ The end -------------------------------
3837 AdvApp2Var_SysBase::mcrdelt_(&c__8, &iszwr, wrkar, &iofwr, &ier);
3842 AdvApp2Var_SysBase::maermsg_("MMA2CE1", iercod, 7L);
3844 AdvApp2Var_SysBase::mgsomsg_("MMA2CE1", 7L);
3849 //=======================================================================
3850 //function : mma2ce2_
3852 //=======================================================================
3853 int mma2ce2_(integer *numdec,
3888 /* System generated locals */
3889 integer sosotb_dim1, sosotb_dim2, sosotb_offset, disotb_dim1, disotb_dim2,
3890 disotb_offset, soditb_dim1, soditb_dim2, soditb_offset,
3891 diditb_dim1, diditb_dim2, diditb_offset, gssutb_dim1, gssvtb_dim1,
3892 chpair_dim1, chpair_dim2, chpair_offset, chimpr_dim1,
3893 chimpr_dim2, chimpr_offset, patjac_dim1, patjac_dim2,
3894 patjac_offset, vecerr_dim1, vecerr_offset, i__1, i__2, i__3, i__4;
3896 /* Local variables */
3897 static logical ldbg;
3898 static integer idim, igsu, minu, minv, maxu, maxv, igsv;
3899 static doublereal vaux[3];
3900 static integer i2rdu, i2rdv, ndses, nd, ii, jj, kk, nu, nv;
3901 static doublereal zu, zv;
3902 static integer nu1, nv1;
3904 /* **********************************************************************
3908 /* Calculation of coefficients of polynomial approximation of degree */
3909 /* (NDJACU,NDJACV) of a function F(u,v), starting from its */
3910 /* discretization on roots of Legendre polynom of degree */
3911 /* NBPNTU by U and NBPNTV by V. */
3915 /* TOUS,AB_SPECIFI::FONCTION&,APPROXIMATION,&COEFFICIENT,&POLYNOME */
3917 /* INPUT ARGUMENTS : */
3918 /* ------------------ */
3919 /* NUMDEC: Indicates if it is POSSIBLE to cut function F(u,v). */
3920 /* = 5, It is POSSIBLE to cut by U or by V or in both directions simultaneously. */
3921 /* = 4, It is POSSIBLE to cut by U or by V BUT NOT in both */
3922 /* directions simultaneously (cutting by V is preferable). */
3923 /* = 3, It is POSSIBLE to cut by U or by V BUT NOT in both */
3924 /* directions simultaneously (cutting by U is preferable). */
3925 /* = 2, It is POSSIBLE to cut only by V (i.e. insert parameter */
3926 /* of cutting Vj). */
3927 /* = 1, It is POSSIBLE to cut only by U (i.e. insert parameter */
3928 /* of cutting Ui). */
3929 /* = 0, It is not POSSIBLE to cut anything */
3930 /* NDIMEN: Total dimension of the space. */
3931 /* NBSESP: Nb of independent sub-spaces on which the errors are calculated. */
3932 /* NDIMSE: Table of dimensions of each of sub-spaces. */
3933 /* NDMINU: Minimum degree by U to be preserved for the approximation. */
3934 /* NDMINV: Minimum degree by V to be preserved for the approximation. */
3935 /* NDGULI: Limit of nb of coefficients by U of the solution. */
3936 /* NDGVLI: Limit of nb of coefficients by V of the solution. */
3937 /* NDJACU: Max degree of the polynom of approximation by U. */
3938 /* The representation in the orthogonal base starts from degree */
3939 /* 0 to degree NDJACU-2*(IORDRU+1). The polynomial base is the base of */
3940 /* Jacobi of order -1 (Legendre), 0, 1 or 2. */
3941 /* It is required that 2*IORDRU+1 <= NDMINU <= NDGULI < NDJACU */
3942 /* NDJACV: Max degree of the polynom of approximation by V. */
3943 /* The representation in the orthogonal base starts from degree */
3944 /* 0 to degree NDJACV-2*(IORDRV+1). The polynomial base is */
3945 /* the base of Jacobi of order -1 (Legendre), 0, 1 or 2 */
3946 /* It is required that 2*IORDRV+1 <= NDMINV <= NDGVLI < NDJACV */
3947 /* IORDRU: Order of the Jacobi base (-1,0,1 or 2) by U. Corresponds */
3948 /* to the step of constraints C0, C1 or C2. */
3949 /* IORDRV: Order of the Jacobi base (-1,0,1 or 2) by U. Corresponds */
3950 /* to the step of constraints C0, C1 or C2. */
3951 /* NBPNTU: Degree of Legendre polynom on the roots which of are */
3952 /* calculated the coefficients of integration by u */
3953 /* by Gauss method. It is required that NBPNTU = 30, 40, */
3954 /* 50 or 61 and NDJACU-2*(IORDRU+1) < NBPNTU. */
3955 /* NBPNTV: Degree of Legendre polynom on the roots which of are */
3956 /* calculated the coefficients of integration by u */
3957 /* by Gauss method. It is required that NBPNTV = 30, 40, */
3958 /* 50 or 61 and NDJACV-2*(IORDRV+1) < NBPNTV. */
3959 /* EPSAPR: Table of NBSESP tolerances imposed on each sub-spaces. */
3960 /* SOSOTB: Table of F(ui,vj) + F(ui,-vj) + F(-ui,vj) + F(-ui,-vj) */
3961 /* with ui and vj - positive roots of the Legendre polynom */
3962 /* of degree NBPNTU and NBPNTV respectively. Additionally, */
3963 /* table SOSOTB(0,j) contains F(0,vj) + F(0,-vj), */
3964 /* table SOSOTB(i,0) contains F(ui,0) + F(-ui,0) and */
3965 /* SOSOTB(0,0) contains F(0,0). */
3966 /* DISOTB: Table of F(ui,vj) + F(ui,-vj) - F(-ui,vj) - F(-ui,-vj) */
3967 /* with ui and vj positive roots of Legendre polynom */
3968 /* of degree NBPNTU and NBPNTV respectively. */
3969 /* SODITB: Table of F(ui,vj) - F(ui,-vj) + F(-ui,vj) - F(-ui,-vj) */
3970 /* with ui and vj positive roots of Legendre polynom */
3971 /* of degree NBPNTU and NBPNTV respectively. */
3972 /* DIDITB: Table of F(ui,vj) - F(ui,-vj) - F(-ui,vj) + F(-ui,-vj) */
3973 /* with ui and vj positive roots of Legendre polynom */
3974 /* of degree NBPNTU and NBPNTV respectively. Additionally, */
3975 /* table DIDITB(0,j) contains F(0,vj) - F(0,-vj), */
3976 /* and table DIDITB(i,0) contains F(ui,0) - F(-ui,0). */
3977 /* GSSUTB: Table of coefficients of integration by Gauss method */
3978 /* by U: i varies from 0 to NBPNTU/2 and k varies from 0 to */
3979 /* NDJACU-2*(IORDRU+1). */
3980 /* GSSVTB: Table of coefficients of integration by Gauss method */
3981 /* by V: i varies from 0 to NBPNTV/2 and k varies from 0 to */
3982 /* NDJACV-2*(IORDRV+1). */
3983 /* XMAXJU: Maximum value of Jacobi polynoms of order IORDRU, */
3984 /* from degree 0 to degree NDJACU - 2*(IORDRU+1) */
3985 /* XMAXJV: Maximum value of Jacobi polynoms of order IORDRV, */
3986 /* from degree 0 to degree NDJACV - 2*(IORDRV+1) */
3988 /* OUTPUT ARGUMENTS : */
3989 /* ------------------- */
3990 /* VECERR: Auxiliary table. */
3991 /* CHPAIR: Auxiliary table of terms connected to degree NDJACU by U */
3992 /* to calculate the coeff. of approximation of EVEN degree by V. */
3993 /* CHIMPR: Auxiliary table of terms connected to degree NDJACU by U */
3994 /* to calculate the coeff. of approximation of UNEVEN degree by V. */
3995 /* PATJAC: Table of coefficients of polynom P(u,v) of approximation */
3996 /* of F(u,v) with eventually taking into account of */
3997 /* constraints. P(u,v) is of degree (NDJACU,NDJACV). */
3998 /* This table contains other coeff if ITYDEC = 0. */
3999 /* ERRMAX: For 1<=i<=NBSESP, ERRMAX(i) contains max errors */
4000 /* on each of sub-spaces SI ITYDEC = 0. */
4001 /* ERRMOY: Contains average errors for each of NBSESP sub-spaces SI ITYDEC = 0. */
4002 /* NDEGPU: Degree by U for square PATJAC. Valable if ITYDEC=0. */
4003 /* NDEGPV: Degree by V for square PATJAC. Valable if ITYDEC=0. */
4004 /* ITYDEC: Shows if it is NECESSARY to cut again function F(u,v). */
4005 /* = 0, it is not NECESSARY to cut anything, PATJAC is OK. */
4006 /* = 1, it is NECESSARY to cut only by U (i.e. insert parameter of cutting Ui). */
4007 /* = 2, it is NECESSARY to cut only by V (i.e. insert parameter of cutting Vj). */
4008 /* = 3, it is NECESSARY to cut both by U AND by V. */
4009 /* IERCOD: Error code. */
4010 /* = 0, Everything is OK. */
4011 /* = -1, There is the best possible solution, but the */
4012 /* user tolerance is not satisfactory (3*only) */
4013 /* = 1, Incoherent entries. */
4015 /* COMMONS USED : */
4016 /* ---------------- */
4018 /* REFERENCES CALLED : */
4019 /* --------------------- */
4021 /* DESCRIPTION/NOTES/LIMITATIONS : */
4023 /* **********************************************************************
4025 /* Name of the routine */
4028 /* --------------------------- Initialisations --------------------------
4031 /* Parameter adjustments */
4032 vecerr_dim1 = *ndimen;
4033 vecerr_offset = vecerr_dim1 + 1;
4034 vecerr -= vecerr_offset;
4039 patjac_dim1 = *ndjacu + 1;
4040 patjac_dim2 = *ndjacv + 1;
4041 patjac_offset = patjac_dim1 * patjac_dim2;
4042 patjac -= patjac_offset;
4043 gssutb_dim1 = *nbpntu / 2 + 1;
4044 chimpr_dim1 = *nbpntv / 2;
4045 chimpr_dim2 = *ndjacu - ((*iordru + 1) << 1) + 1;
4046 chimpr_offset = chimpr_dim1 * chimpr_dim2 + 1;
4047 chimpr -= chimpr_offset;
4048 chpair_dim1 = *nbpntv / 2 + 1;
4049 chpair_dim2 = *ndjacu - ((*iordru + 1) << 1) + 1;
4050 chpair_offset = chpair_dim1 * chpair_dim2;
4051 chpair -= chpair_offset;
4052 gssvtb_dim1 = *nbpntv / 2 + 1;
4053 diditb_dim1 = *nbpntu / 2 + 1;
4054 diditb_dim2 = *nbpntv / 2 + 1;
4055 diditb_offset = diditb_dim1 * diditb_dim2;
4056 diditb -= diditb_offset;
4057 soditb_dim1 = *nbpntu / 2;
4058 soditb_dim2 = *nbpntv / 2;
4059 soditb_offset = soditb_dim1 * (soditb_dim2 + 1) + 1;
4060 soditb -= soditb_offset;
4061 disotb_dim1 = *nbpntu / 2;
4062 disotb_dim2 = *nbpntv / 2;
4063 disotb_offset = disotb_dim1 * (disotb_dim2 + 1) + 1;
4064 disotb -= disotb_offset;
4065 sosotb_dim1 = *nbpntu / 2 + 1;
4066 sosotb_dim2 = *nbpntv / 2 + 1;
4067 sosotb_offset = sosotb_dim1 * sosotb_dim2;
4068 sosotb -= sosotb_offset;
4071 ldbg = AdvApp2Var_SysBase::mnfndeb_() >= 3;
4073 AdvApp2Var_SysBase::mgenmsg_("MMA2CE2", 7L);
4075 /* --> A priori everything is OK */
4077 /* --> test of inputs */
4078 if (*numdec < 0 || *numdec > 5) {
4081 if ((*iordru << 1) + 1 > *ndminu) {
4084 if (*ndminu > *ndguli) {
4087 if (*ndguli >= *ndjacu) {
4090 if ((*iordrv << 1) + 1 > *ndminv) {
4093 if (*ndminv > *ndgvli) {
4096 if (*ndgvli >= *ndjacv) {
4099 /* --> A priori, no cuts to be done */
4101 /* --> Min. degrees to return: NDMINU,NDMINV */
4104 /* --> For the moment, max errors are null */
4105 AdvApp2Var_SysBase::mvriraz_(nbsesp, (char *)&errmax[1]);
4107 AdvApp2Var_SysBase::mvriraz_(&nd, (char *)&vecerr[vecerr_offset]);
4108 /* --> and the square, too. */
4109 nd = (*ndjacu + 1) * (*ndjacv + 1) * *ndimen;
4110 AdvApp2Var_SysBase::mvriraz_(&nd, (char *)&patjac[patjac_offset]);
4112 i2rdu = (*iordru + 1) << 1;
4113 i2rdv = (*iordrv + 1) << 1;
4115 /* **********************************************************************
4117 /* -------------------- HERE IT IS POSSIBLE TO CUT ----------------------
4119 /* **********************************************************************
4122 if (*numdec > 0 && *numdec <= 5) {
4124 /* ******************************************************************
4126 /* ---------------------- Calculate coeff of zone 4 -------------
4140 /* ---------------- Calculate the terms connected to degree by U ---------
4144 for (nd = 1; nd <= i__1; ++nd) {
4146 for (kk = minu; kk <= i__2; ++kk) {
4148 mma2cfu_(&kk, nbpntu, nbpntv, &sosotb[nd * sosotb_dim2 *
4149 sosotb_dim1], &disotb[(nd * disotb_dim2 + 1) *
4150 disotb_dim1 + 1], &soditb[(nd * soditb_dim2 + 1) *
4151 soditb_dim1 + 1], &diditb[nd * diditb_dim2 *
4152 diditb_dim1], &gssutb[igsu * gssutb_dim1], &chpair[(
4153 igsu + nd * chpair_dim2) * chpair_dim1], &chimpr[(
4154 igsu + nd * chimpr_dim2) * chimpr_dim1 + 1]);
4160 /* ------------------- Calculate the coefficients of PATJAC ------------
4163 igsu = minu - i2rdu;
4165 for (jj = minv; jj <= i__1; ++jj) {
4168 for (nd = 1; nd <= i__2; ++nd) {
4169 mma2cfv_(&jj, &minu, &maxu, nbpntv, &gssvtb[igsv *
4170 gssvtb_dim1], &chpair[(igsu + nd * chpair_dim2) *
4171 chpair_dim1], &chimpr[(igsu + nd * chimpr_dim2) *
4172 chimpr_dim1 + 1], &patjac[minu + (jj + nd *
4173 patjac_dim2) * patjac_dim1]);
4177 /* ----- Contribution of calculated terms to the approximation error */
4178 /* for terms (I,J) with MINU <= I <= MAXU, J fixe. */
4182 for (nd = 1; nd <= i__2; ++nd) {
4184 mma2er1_(ndjacu, ndjacv, &ndses, &minu, &maxu, &jj, &jj,
4185 iordru, iordrv, xmaxju, xmaxjv, &patjac[idim *
4186 patjac_dim2 * patjac_dim1], &vecerr[vecerr_dim1 + 1],
4187 &vecerr[nd + (vecerr_dim1 << 2)]);
4188 if (vecerr[nd + (vecerr_dim1 << 2)] > epsapr[nd]) {
4197 /* ******************************************************************
4199 /* ---------------------- Calculate the coeff of zone 2 -------------
4202 minu = (*iordru + 1) << 1;
4207 /* --> If zone 2 is empty, pass to zone 3. */
4208 /* VECERR(ND,2) was already set to zero. */
4213 /* ---------------- Calculate the terms connected to degree by U ------------
4217 for (nd = 1; nd <= i__1; ++nd) {
4219 for (kk = minu; kk <= i__2; ++kk) {
4221 mma2cfu_(&kk, nbpntu, nbpntv, &sosotb[nd * sosotb_dim2 *
4222 sosotb_dim1], &disotb[(nd * disotb_dim2 + 1) *
4223 disotb_dim1 + 1], &soditb[(nd * soditb_dim2 + 1) *
4224 soditb_dim1 + 1], &diditb[nd * diditb_dim2 *
4225 diditb_dim1], &gssutb[igsu * gssutb_dim1], &chpair[(
4226 igsu + nd * chpair_dim2) * chpair_dim1], &chimpr[(
4227 igsu + nd * chimpr_dim2) * chimpr_dim1 + 1]);
4233 /* ------------------- Calculate the coefficients of PATJAC ------------
4236 igsu = minu - i2rdu;
4238 for (jj = minv; jj <= i__1; ++jj) {
4241 for (nd = 1; nd <= i__2; ++nd) {
4242 mma2cfv_(&jj, &minu, &maxu, nbpntv, &gssvtb[igsv *
4243 gssvtb_dim1], &chpair[(igsu + nd * chpair_dim2) *
4244 chpair_dim1], &chimpr[(igsu + nd * chimpr_dim2) *
4245 chimpr_dim1 + 1], &patjac[minu + (jj + nd *
4246 patjac_dim2) * patjac_dim1]);
4252 /* -----Contribution of calculated terms to the approximation error */
4253 /* for terms (I,J) with MINU <= I <= MAXU, MINV <= J <= MAXV */
4257 for (nd = 1; nd <= i__1; ++nd) {
4259 mma2er1_(ndjacu, ndjacv, &ndses, &minu, &maxu, &minv, &maxv,
4260 iordru, iordrv, xmaxju, xmaxjv, &patjac[idim *
4261 patjac_dim2 * patjac_dim1], &vecerr[vecerr_dim1 + 1], &
4262 vecerr[nd + (vecerr_dim1 << 1)]);
4267 /* ******************************************************************
4269 /* ---------------------- Calculation of coeff of zone 3 -------------
4275 minv = (*iordrv + 1) << 1;
4278 /* -> If zone 3 is empty, pass to the test of cutting. */
4279 /* VECERR(ND,3) was already set to zero */
4284 /* ----------- The terms connected to the degree by U are already calculated -----
4286 /* ------------------- Calculation of coefficients of PATJAC ------------
4289 igsu = minu - i2rdu;
4291 for (jj = minv; jj <= i__1; ++jj) {
4294 for (nd = 1; nd <= i__2; ++nd) {
4295 mma2cfv_(&jj, &minu, &maxu, nbpntv, &gssvtb[igsv *
4296 gssvtb_dim1], &chpair[(igsu + nd * chpair_dim2) *
4297 chpair_dim1], &chimpr[(igsu + nd * chimpr_dim2) *
4298 chimpr_dim1 + 1], &patjac[minu + (jj + nd *
4299 patjac_dim2) * patjac_dim1]);
4305 /* ----- Contribution of calculated terms to the approximation error
4306 /* for terms (I,J) with MINU <= I <= MAXU, MINV <= J <= MAXV. */
4310 for (nd = 1; nd <= i__1; ++nd) {
4312 mma2er1_(ndjacu, ndjacv, &ndses, &minu, &maxu, &minv, &maxv,
4313 iordru, iordrv, xmaxju, xmaxjv, &patjac[idim *
4314 patjac_dim2 * patjac_dim1], &vecerr[vecerr_dim1 + 1], &
4315 vecerr[nd + vecerr_dim1 * 3]);
4320 /* ******************************************************************
4322 /* --------------------------- Tests of cutting ---------------------
4327 for (nd = 1; nd <= i__1; ++nd) {
4328 vaux[0] = vecerr[nd + (vecerr_dim1 << 1)];
4329 vaux[1] = vecerr[nd + (vecerr_dim1 << 2)];
4330 vaux[2] = vecerr[nd + vecerr_dim1 * 3];
4332 errmax[nd] = AdvApp2Var_MathBase::mzsnorm_(&ii, vaux);
4333 if (errmax[nd] > epsapr[nd]) {
4335 zv = AdvApp2Var_MathBase::mzsnorm_(&ii, vaux);
4336 zu = AdvApp2Var_MathBase::mzsnorm_(&ii, &vaux[1]);
4337 if (zu > epsapr[nd] && zv > epsapr[nd]) {
4349 /* ******************************************************************
4351 /* --- OK, the square is valid, the coeff of zone 1 are calculated
4354 minu = (*iordru + 1) << 1;
4356 minv = (*iordrv + 1) << 1;
4359 /* --> If zone 1 is empty, pass to the calculation of Max and Average error. */
4360 if (minu > maxu || minv > maxv) {
4364 /* ----------- The terms connected to degree by U are already calculated -----
4366 /* ------------------- Calculate the coefficients of PATJAC ------------
4369 igsu = minu - i2rdu;
4371 for (jj = minv; jj <= i__1; ++jj) {
4374 for (nd = 1; nd <= i__2; ++nd) {
4375 mma2cfv_(&jj, &minu, &maxu, nbpntv, &gssvtb[igsv *
4376 gssvtb_dim1], &chpair[(igsu + nd * chpair_dim2) *
4377 chpair_dim1], &chimpr[(igsu + nd * chimpr_dim2) *
4378 chimpr_dim1 + 1], &patjac[minu + (jj + nd *
4379 patjac_dim2) * patjac_dim1]);
4385 /* --------------- Now the degree is maximally lowered --------
4390 i__1 = 1, i__2 = (*iordru << 1) + 1, i__1 = advapp_max(i__1,i__2);
4391 minu = advapp_max(i__1,*ndminu);
4394 i__1 = 1, i__2 = (*iordrv << 1) + 1, i__1 = advapp_max(i__1,i__2);
4395 minv = advapp_max(i__1,*ndminv);
4399 for (nd = 1; nd <= i__1; ++nd) {
4401 if (maxu >= (*iordru + 1) << 1 && maxv >= (*iordrv + 1) << 1) {
4402 mma2er2_(ndjacu, ndjacv, &ndses, &minu, &maxu, &minv, &maxv,
4403 iordru, iordrv, xmaxju, xmaxjv, &patjac[idim *
4404 patjac_dim2 * patjac_dim1], &epsapr[nd], &vecerr[
4405 vecerr_dim1 + 1], &errmax[nd], &nu, &nv);
4413 /* --> Calculate the average error. */
4414 mma2moy_(ndjacu, ndjacv, &ndses, &nu1, ndjacu, &nv1, ndjacv,
4415 iordru, iordrv, &patjac[idim * patjac_dim2 * patjac_dim1],
4418 /* --> Set to 0.D0 the rejected coeffs. */
4419 i__2 = idim + ndses - 1;
4420 for (ii = idim; ii <= i__2; ++ii) {
4422 for (jj = nv1; jj <= i__3; ++jj) {
4424 for (kk = nu1; kk <= i__4; ++kk) {
4425 patjac[kk + (jj + ii * patjac_dim2) * patjac_dim1] =
4434 /* --> Return the nb of coeffs of approximation. */
4435 *ndegpu = advapp_max(*ndegpu,nu);
4436 *ndegpv = advapp_max(*ndegpv,nv);
4441 /* ******************************************************************
4443 /* -------------------- IT IS NOT POSSIBLE TO CUT -------------------
4445 /* ******************************************************************
4449 minu = (*iordru + 1) << 1;
4451 minv = (*iordrv + 1) << 1;
4454 /* ---------------- Calculate the terms connected to the degree by U ------------
4458 for (nd = 1; nd <= i__1; ++nd) {
4460 for (kk = minu; kk <= i__2; ++kk) {
4462 mma2cfu_(&kk, nbpntu, nbpntv, &sosotb[nd * sosotb_dim2 *
4463 sosotb_dim1], &disotb[(nd * disotb_dim2 + 1) *
4464 disotb_dim1 + 1], &soditb[(nd * soditb_dim2 + 1) *
4465 soditb_dim1 + 1], &diditb[nd * diditb_dim2 *
4466 diditb_dim1], &gssutb[igsu * gssutb_dim1], &chpair[(
4467 igsu + nd * chpair_dim2) * chpair_dim1], &chimpr[(
4468 igsu + nd * chimpr_dim2) * chimpr_dim1 + 1]);
4472 /* ---------------------- Calculate all coefficients -------
4475 igsu = minu - i2rdu;
4477 for (jj = minv; jj <= i__2; ++jj) {
4479 mma2cfv_(&jj, &minu, &maxu, nbpntv, &gssvtb[igsv *
4480 gssvtb_dim1], &chpair[(igsu + nd * chpair_dim2) *
4481 chpair_dim1], &chimpr[(igsu + nd * chimpr_dim2) *
4482 chimpr_dim1 + 1], &patjac[minu + (jj + nd *
4483 patjac_dim2) * patjac_dim1]);
4489 /* ----- Contribution of calculated terms to the approximation error
4490 /* for terms (I,J) with MINU <= I <= MAXU, MINV <= J <= MAXV */
4494 for (nd = 1; nd <= i__1; ++nd) {
4496 minu = (*iordru + 1) << 1;
4500 mma2er1_(ndjacu, ndjacv, &ndses, &minu, &maxu, &minv, &maxv,
4501 iordru, iordrv, xmaxju, xmaxjv, &patjac[idim *
4502 patjac_dim2 * patjac_dim1], &vecerr[vecerr_dim1 + 1], &
4506 minv = (*iordrv + 1) << 1;
4509 mma2er1_(ndjacu, ndjacv, &ndses, &minu, &maxu, &minv, &maxv,
4510 iordru, iordrv, xmaxju, xmaxjv, &patjac[idim *
4511 patjac_dim2 * patjac_dim1], &vecerr[vecerr_dim1 + 1],
4515 /* ---------------------------- IF ERRMAX > EPSAPR, stop --------
4518 if (errmax[nd] > epsapr[nd]) {
4523 /* ------------- Otherwise, try to remove again the coeff
4528 i__2 = 1, i__3 = (*iordru << 1) + 1, i__2 = advapp_max(i__2,i__3);
4529 minu = advapp_max(i__2,*ndminu);
4532 i__2 = 1, i__3 = (*iordrv << 1) + 1, i__2 = advapp_max(i__2,i__3);
4533 minv = advapp_max(i__2,*ndminv);
4535 if (maxu >= (*iordru + 1) << 1 && maxv >= (*iordrv + 1) << 1) {
4536 mma2er2_(ndjacu, ndjacv, &ndses, &minu, &maxu, &minv, &
4537 maxv, iordru, iordrv, xmaxju, xmaxjv, &patjac[
4538 idim * patjac_dim2 * patjac_dim1], &epsapr[nd], &
4539 vecerr[vecerr_dim1 + 1], &errmax[nd], &nu, &nv);
4546 /* --------------------- Calculate the average error -------------
4551 mma2moy_(ndjacu, ndjacv, &ndses, &nu1, ndjacu, &nv1, ndjacv,
4552 iordru, iordrv, &patjac[idim * patjac_dim2 * patjac_dim1],
4555 /* --------------------- Set to 0.D0 the rejected coeffs ----------
4558 i__2 = idim + ndses - 1;
4559 for (ii = idim; ii <= i__2; ++ii) {
4561 for (jj = nv1; jj <= i__3; ++jj) {
4563 for (kk = nu1; kk <= i__4; ++kk) {
4564 patjac[kk + (jj + ii * patjac_dim2) * patjac_dim1] =
4573 /* --------------- Return the nb of coeff of approximation ---
4576 *ndegpu = advapp_max(*ndegpu,nu);
4577 *ndegpv = advapp_max(*ndegpv,nv);
4585 /* ------------------------------ The end -------------------------------
4587 /* --> Error in inputs */
4592 /* --------- Management of cuts, it is required 0 < NUMDEC <= 5 -------
4595 /* --> Here it is possible and necessary to cut, choose by U if it is possible */
4597 if (*numdec <= 0 || *numdec > 5) {
4606 /* --> Here it is possible and necessary to cut, choose by U if it is possible */
4608 if (*numdec <= 0 || *numdec > 5) {
4617 /* --> Here it is possible and necessary to cut, choose by 4 if it is possible */
4619 if (*numdec <= 0 || *numdec > 5) {
4624 } else if (*numdec == 2 || *numdec == 4) {
4626 } else if (*numdec == 1 || *numdec == 3) {
4634 AdvApp2Var_SysBase::maermsg_("MMA2CE2", iercod, 7L);
4636 AdvApp2Var_SysBase::mgsomsg_("MMA2CE2", 7L);
4641 //=======================================================================
4642 //function : mma2cfu_
4644 //=======================================================================
4645 int mma2cfu_(integer *ndujac,
4657 /* System generated locals */
4658 integer sosotb_dim1, disotb_dim1, disotb_offset, soditb_dim1,
4659 soditb_offset, diditb_dim1, i__1, i__2;
4661 /* Local variables */
4662 static logical ldbg;
4663 static integer nptu2, nptv2, ii, jj;
4664 static doublereal bid0, bid1, bid2;
4666 /* **********************************************************************
4671 /* Calculate the terms connected to degree NDUJAC by U of the polynomial approximation */
4672 /* of function F(u,v), starting from its discretisation
4673 /* on the roots of Legendre polynom of degree */
4674 /* NBPNTU by U and NBPNTV by V. */
4678 /* FONCTION,APPROXIMATION,COEFFICIENT,POLYNOME */
4680 /* INPUT ARGUMENTSE : */
4681 /* ------------------ */
4682 /* NDUJAC: Fixed degree by U for which the terms */
4683 /* allowing to obtain the Legendre or Jacobi coeff*/
4684 /* of even or uneven degree by V are calculated. */
4685 /* NBPNTU: Degree of Legendre polynom on the roots which of */
4686 /* the coefficients of integration by U are calculated */
4687 /* by Gauss method. It is required that NBPNTU = 30, 40, 50 or 61. */
4688 /* NBPNTV: Degree of Legendre polynom on the roots which of */
4689 /* the coefficients of integration by V are calculated */
4690 /* by Gauss method. It is required that NBPNTV = 30, 40, 50 or 61. */
4691 /* SOSOTB: Table of F(ui,vj) + F(ui,-vj) + F(-ui,vj) + F(-ui,-vj) */
4692 /* with ui and vj positive roots of Legendre polynom */
4693 /* of degree NBPNTU and NBPNTV respectively. Moreover, */
4694 /* table SOSOTB(0,j) contains F(0,vj) + F(0,-vj), */
4695 /* table SOSOTB(i,0) contains F(ui,0) + F(-ui,0) and */
4696 /* SOSOTB(0,0) contains F(0,0). */
4697 /* DISOTB: Table of F(ui,vj) + F(ui,-vj) - F(-ui,vj) - F(-ui,-vj) */
4698 /* with ui and vj positive roots of Legendre polynom */
4699 /* of degree NBPNTU and NBPNTV respectively. */
4700 /* SODITB: Table of F(ui,vj) - F(ui,-vj) + F(-ui,vj) - F(-ui,-vj) */
4701 /* with ui and vj positive roots of Legendre polynom */
4702 /* of degree NBPNTU and NBPNTV respectively. */
4703 /* DIDITB: Table of F(ui,vj) - F(ui,-vj) - F(-ui,vj) + F(-ui,-vj) */
4704 /* avec ui and vj positive roots of Legendre polynom */
4705 /* of degree NBPNTU and NBPNTV respectively. Moreover, */
4706 /* table DIDITB(0,j) contains F(0,vj) - F(0,-vj), */
4707 /* and table DIDITB(i,0) contains F(ui,0) - F(-ui,0). */
4708 /* GSSUTB: Table of coefficients of integration by Gauss method */
4709 /* Gauss by U for fixed NDUJAC : i varies from 0 to NBPNTU/2. */
4711 /* OUTPUT ARGUMENTS : */
4712 /* ------------------- */
4713 /* CHPAIR: Table of terms connected to degree NDUJAC by U to calculate the */
4714 /* coeff. of the approximation of EVEN degree by V. */
4715 /* CHIMPR: Table of terms connected to degree NDUJAC by U to calculate */
4716 /* the coeff. of approximation of UNEVEN degree by V. */
4718 /* COMMONS USED : */
4719 /* ---------------- */
4721 /* REFERENCES CALLED : */
4722 /* ----------------------- */
4724 /* DESCRIPTION/NOTES/LIMITATIONS : */
4725 /* ----------------------------------- */
4729 /* **********************************************************************
4731 /* Name of the routine */
4734 /* --------------------------- Initialisations --------------------------
4737 /* Parameter adjustments */
4739 diditb_dim1 = *nbpntu / 2 + 1;
4740 soditb_dim1 = *nbpntu / 2;
4741 soditb_offset = soditb_dim1 + 1;
4742 soditb -= soditb_offset;
4743 disotb_dim1 = *nbpntu / 2;
4744 disotb_offset = disotb_dim1 + 1;
4745 disotb -= disotb_offset;
4746 sosotb_dim1 = *nbpntu / 2 + 1;
4749 ldbg = AdvApp2Var_SysBase::mnfndeb_() >= 3;
4751 AdvApp2Var_SysBase::mgenmsg_("MMA2CFU", 7L);
4754 nptu2 = *nbpntu / 2;
4755 nptv2 = *nbpntv / 2;
4757 /* **********************************************************************
4759 /* CALCULATE COEFFICIENTS BY U */
4761 /* ----------------- Calculate coefficients of even degree --------------
4764 if (*ndujac % 2 == 0) {
4766 for (jj = 1; jj <= i__1; ++jj) {
4770 for (ii = 1; ii <= i__2; ++ii) {
4772 bid1 += sosotb[ii + jj * sosotb_dim1] * bid0;
4773 bid2 += soditb[ii + jj * soditb_dim1] * bid0;
4781 /* --------------- Calculate coefficients of uneven degree ----------
4786 for (jj = 1; jj <= i__1; ++jj) {
4790 for (ii = 1; ii <= i__2; ++ii) {
4792 bid1 += disotb[ii + jj * disotb_dim1] * bid0;
4793 bid2 += diditb[ii + jj * diditb_dim1] * bid0;
4802 /* ------- Add terms connected to the supplementary root (0.D0) ------
4803 /* ----------- of Legendre polynom of uneven degree NBPNTU -----------
4805 /* --> Only even NDUJAC terms are modified as GSSUTB(0) = 0 */
4806 /* when NDUJAC is uneven. */
4808 if (*nbpntu % 2 != 0 && *ndujac % 2 == 0) {
4811 for (jj = 1; jj <= i__1; ++jj) {
4812 chpair[jj] += sosotb[jj * sosotb_dim1] * bid0;
4813 chimpr[jj] += diditb[jj * diditb_dim1] * bid0;
4818 /* ------ Calculate the terms connected to supplementary roots (0.D0) ------
4820 /* ----------- of Legendre polynom of uneven degree NBPNTV -----------
4823 if (*nbpntv % 2 != 0) {
4824 /* --> Only CHPAIR terms are calculated as GSSVTB(0,IH-IDEBV)=0
4826 /* when IH is uneven (see MMA2CFV). */
4828 if (*ndujac % 2 == 0) {
4831 for (ii = 1; ii <= i__1; ++ii) {
4832 bid1 += sosotb[ii] * gssutb[ii];
4839 for (ii = 1; ii <= i__1; ++ii) {
4840 bid1 += diditb[ii] * gssutb[ii];
4845 if (*nbpntu % 2 != 0) {
4846 chpair[0] += sosotb[0] * gssutb[0];
4850 /* ------------------------------ The end -------------------------------
4854 AdvApp2Var_SysBase::mgsomsg_("MMA2CFU", 7L);
4859 //=======================================================================
4860 //function : mma2cfv_
4862 //=======================================================================
4863 int mma2cfv_(integer *ndvjac,
4873 /* System generated locals */
4874 integer chpair_dim1, chpair_offset, chimpr_dim1, chimpr_offset,
4875 patjac_offset, i__1, i__2;
4877 /* Local variables */
4878 static logical ldbg;
4879 static integer nptv2, ii, jj;
4880 static doublereal bid1;
4882 /* **********************************************************************
4887 /* Calculate the coefficients of polynomial approximation of F(u,v)
4888 /* of degree NDVJAC by V and of degree by U varying from MINDGU to MAXDGU.
4893 /* FONCTION,APPROXIMATION,COEFFICIENT,POLYNOME */
4895 /* INPUT ARGUMENTS : */
4896 /* ------------------ */
4898 /* NDVJAC: Degree of the polynom of approximation by V. */
4899 /* The representation in the orthogonal base starts from degre 0.
4900 /* The polynomial base is the base of Jacobi of order -1 */
4901 /* (Legendre), 0, 1 or 2 */
4902 /* MINDGU: Degree minimum by U of coeff. to calculate. */
4903 /* MAXDGU: Degree maximum by U of coeff. to calculate. */
4904 /* NBPNTV: Degree of the Legendre polynom on the roots which of */
4905 /* the coefficients of integration by V are calculated */
4906 /* by Gauss method. It is reqired that NBPNTV = 30, 40, 50 or 61 and NDVJAC < NBPNTV. */
4907 /* GSSVTB: Table of coefficients of integration by Gauss method */
4908 /* by V for NDVJAC fixed: j varies from 0 to NBPNTV/2. */
4909 /* CHPAIR: Table of terms connected to degrees from MINDGU to MAXDGU by U to
4910 /* calculate the coeff. of approximation of EVEN degree NDVJAC by V. */
4911 /* CHIMPR: Table of terms connected to degrees from MINDGU to MAXDGU by U to
4912 /* calculate the coeff. of approximation of UNEVEN degree NDVJAC by V. */
4914 /* OUTPUT ARGUMENTS : */
4915 /* ------------------- */
4916 /* PATJAC: Table of coefficients by U of the polynom of approximation */
4917 /* P(u,v) of degree MINDGU to MAXDGU by U and NDVJAC by V. */
4919 /* COMMONS USED : */
4920 /* -------------- */
4922 /* REFERENCES CALLED : */
4923 /* --------------------- */
4925 /* DESCRIPTION/NOTES/LIMITATIONS : */
4926 /* ------------------------------- */
4928 /* **********************************************************************
4930 /* Name of the routine */
4933 /* --------------------------- Initialisations --------------------------
4936 /* Parameter adjustments */
4937 patjac_offset = *mindgu;
4938 patjac -= patjac_offset;
4939 chimpr_dim1 = *nbpntv / 2;
4940 chimpr_offset = chimpr_dim1 * *mindgu + 1;
4941 chimpr -= chimpr_offset;
4942 chpair_dim1 = *nbpntv / 2 + 1;
4943 chpair_offset = chpair_dim1 * *mindgu;
4944 chpair -= chpair_offset;
4947 ldbg = AdvApp2Var_SysBase::mnfndeb_() >= 3;
4949 AdvApp2Var_SysBase::mgenmsg_("MMA2CFV", 7L);
4951 nptv2 = *nbpntv / 2;
4953 /* --------- Calculate the coefficients for even degree NDVJAC ----------
4956 if (*ndvjac % 2 == 0) {
4958 for (ii = *mindgu; ii <= i__1; ++ii) {
4961 for (jj = 1; jj <= i__2; ++jj) {
4962 bid1 += chpair[jj + ii * chpair_dim1] * gssvtb[jj];
4969 /* -------- Calculate the coefficients for uneven degree NDVJAC -----
4974 for (ii = *mindgu; ii <= i__1; ++ii) {
4977 for (jj = 1; jj <= i__2; ++jj) {
4978 bid1 += chimpr[jj + ii * chimpr_dim1] * gssvtb[jj];
4986 /* ------- Add terms connected to the supplementary root (0.D0) ----- */
4987 /* --------of the Legendre polynom of uneven degree NBPNTV --------- */
4989 if (*nbpntv % 2 != 0 && *ndvjac % 2 == 0) {
4992 for (ii = *mindgu; ii <= i__1; ++ii) {
4993 patjac[ii] += bid1 * chpair[ii * chpair_dim1];
4998 /* ------------------------------ The end -------------------------------
5002 AdvApp2Var_SysBase::mgsomsg_("MMA2CFV", 7L);
5007 //=======================================================================
5008 //function : mma2ds1_
5010 //=======================================================================
5011 int AdvApp2Var_ApproxF2var::mma2ds1_(integer *ndimen,
5014 const AdvApp2Var_EvaluatorFunc2Var& foncnp,
5029 /* System generated locals */
5030 integer sosotb_dim1, sosotb_dim2, sosotb_offset, disotb_dim1, disotb_dim2,
5031 disotb_offset, soditb_dim1, soditb_dim2, soditb_offset,
5032 diditb_dim1, diditb_dim2, diditb_offset, fpntab_dim1,
5033 fpntab_offset, i__1;
5035 /* Local variables */
5036 static logical ldbg;
5037 static integer ibid1, ibid2, iuouv, nd;
5038 static integer isz1, isz2;
5040 /* **********************************************************************
5045 /* Discretisation of function F(u,v) on the roots of Legendre polynoms. */
5049 /* FONCTION&,DISCRETISATION,&POINT */
5051 /* INPUT ARGUMENTS : */
5052 /* ------------------ */
5053 /* NDIMEN: Dimension of the space. */
5054 /* UINTFN: Limits of the interval of definition by u of the function */
5055 /* to be processed: (UINTFN(1),UINTFN(2)). */
5056 /* VINTFN: Limits of the interval of definition by v of the function */
5057 /* to be processed: (VINTFN(1),VINTFN(2)). */
5058 /* FONCNP: The NAME of the non-polynomial function to be processed. */
5059 /* NBPNTU: The degree of Legendre polynom on the roots which of */
5060 /* FONCNP is discretized by u. */
5061 /* NBPNTV: The degree of Legendre polynom on the roots which of */
5062 /* FONCNP is discretized by v. */
5063 /* UROOTB: Table of STRICTLY POSITIVE roots of the polynom */
5064 /* of Legendre of degree NBPNTU defined on (-1,1). */
5065 /* VROOTB: Table of STRICTLY POSITIVE roots of the polynom */
5066 /* of Legendre of degree NBPNTV defined on (-1,1). */
5067 /* ISOFAV: Shows the type of iso of F(u,v) to be extracted to improve */
5068 /* the rapidity of calculation (has no influence on the form */
5070 /* = 1, shows that it is necessary to calculate the points of F(u,v) */
5071 /* with fixed u (with NBPNTV values different from v). */
5072 /* = 2, shows that it is necessaty to calculate the points of F(u,v) */
5073 /* with fixed v (with NBPNTU values different from u). */
5074 /* SOSOTB: Preinitialized table (input/output argument). */
5075 /* DISOTB: Preinitialized table (input/output argument). */
5076 /* SODITB: Preinitialized table (input/output argument). */
5077 /* DIDITB: Preinitialized table (input/output argument). */
5079 /* OUTPUT ARGUMENTS : */
5080 /* ------------------- */
5081 /* SOSOTB: Table where the terms */
5082 /* F(ui,vj) + F(ui,-vj) + F(-ui,vj) + F(-ui,-vj) */
5083 /* are added with ui and vj positive roots of Legendre polynom */
5084 /* of degree NBPNTU and NBPNTV respectively. */
5085 /* DISOTB: Table where the terms */
5086 /* F(ui,vj) + F(ui,-vj) - F(-ui,vj) - F(-ui,-vj) */
5087 /* are added with ui and vj positive roots of Legendre polynom */
5088 /* of degree NBPNTU and NBPNTV respectively. */
5089 /* SODITB: Table where the terms */
5090 /* F(ui,vj) - F(ui,-vj) + F(-ui,vj) - F(-ui,-vj) */
5091 /* are added with ui and vj positive roots of Legendre polynom */
5092 /* of degree NBPNTU and NBPNTV respectively. */
5093 /* DIDITB: Table where the terms */
5094 /* F(ui,vj) - F(ui,-vj) - F(-ui,vj) + F(-ui,-vj) */
5095 /* are added with ui and vj positive roots of Legendre polynom */
5096 /* of degree NBPNTU and NBPNTV respectively. */
5097 /* FPNTAB: Auxiliary table. */
5098 /* TTABLE: Auxiliary table. */
5099 /* IERCOD: Error code >100 Pb in the evaluation of FONCNP, */
5100 /* the returned error code is equal to error code of FONCNP + 100. */
5102 /* COMMONS USED : */
5103 /* ---------------- */
5105 /* REFERENCES CALLED : */
5106 /* --------------------- */
5108 /* DESCRIPTION/NOTES/LIMITATIONS : */
5109 /* ----------------------------------- */
5110 /* --> The external function created by the caller of MA2F1K, MA2FDK */
5111 /* where MA2FXK should be in the following form : */
5112 /* SUBROUTINE FONCNP(NDIMEN,UINTFN,VINTFN,ISOFAV,TCONST,NBPTAB */
5113 /* ,TTABLE,IDERIU,IDERIV,PPNTAB,IERCOD) */
5114 /* with the following input arguments : */
5115 /* - NDIMEN is integer defined as the sum of dimensions of */
5116 /* sub-spaces (i.e. total dimension of the problem). */
5117 /* - UINTFN(2) is a table of 2 reals containing the interval */
5118 /* by u where the function to be approximated is defined */
5119 /* (so it is equal to UIFONC). */
5120 /* - VINTFN(2) is a table of 2 reals containing the interval */
5121 /* by v where the function to be approximated is defined */
5122 /* (so it is equal to VIFONC). */
5123 /* - ISOFAV, is 1 if it is necessary to calculate points with constant u, */
5124 /* is 2 if it is necessary to calculate points with constant v. */
5125 /* Any other value is an error. */
5126 /* - TCONST, real, value of the fixed parameter. Takes values */
5127 /* in (UIFONC(1),UIFONC(2)) if ISOFAV = 1 or */
5128 /* ins (VIFONC(1),VIFONC(2)) if ISOFAV = 2. */
5129 /* - NBPTAB, integer. Shows the number of points to be calculated. */
5130 /* - TTABLE, a table of reals NBPTAB. These are the values of */
5131 /* 'free' parameter of discretization (v if IISOFAV=1, */
5132 /* u if IISOFAV=2). */
5133 /* - IDERIU, integer, takes values between 0 (position) */
5134 /* and IORDRE(1) (partial derivative of the function by u */
5135 /* of order IORDRE(1) if IORDRE(1) > 0). */
5136 /* - IDERIV, integer, takes values between 0 (position) */
5137 /* and IORDRE(2) (partial derivative of the function by v */
5138 /* of order IORDRE(2) if IORDRE(2) > 0). */
5139 /* If IDERIU=i and IDERIV=j, FONCNP should calculate the */
5140 /* points of the derivative : */
5147 /* and the output arguments aret : */
5148 /* - FPNTAB(NDIMEN,NBPTAB) contains, at output, the table of */
5149 /* NBPTAB points calculated in FONCNP. */
5150 /* - IERCOD is, at output the error code of FONCNP. This code */
5151 /* (integer) should be strictly positive if there is a problem. */
5153 /* The input arguments SHOULD NOT be modified under FONCNP.
5156 /* -->As FONCNP is not forcedly defined in (-1,1)*(-1,1), the */
5157 /* values of UROOTB and VROOTB are consequently modified. */
5159 /* -->The results of discretisation are ranked in 4 tables */
5160 /* SOSOTB, DISOTB, SODITB and DIDITB to earn time */
5161 /* during the calculation of coefficients of the polynom of approximation. */
5163 /* When NBPNTU is uneven : */
5164 /* table SOSOTB(0,j) contains F(0,vj) + F(0,-vj), */
5165 /* table DIDITB(0,j) contains F(0,vj) - F(0,-vj), */
5166 /* When NBPNTV is uneven : */
5167 /* table SOSOTB(i,0) contains F(ui,0) + F(-ui,0), */
5168 /* table DIDITB(i,0) contains F(ui,0) - F(-ui,0), */
5169 /* When NBPNTU and NBPNTV are uneven : */
5170 /* term SOSOTB(0,0) contains F(0,0). */
5173 /* **********************************************************************
5175 /* Name of the routine */
5178 /* --------------------------- Initialization --------------------------
5181 /* Parameter adjustments */
5182 fpntab_dim1 = *ndimen;
5183 fpntab_offset = fpntab_dim1 + 1;
5184 fpntab -= fpntab_offset;
5188 diditb_dim1 = *nbpntu / 2 + 1;
5189 diditb_dim2 = *nbpntv / 2 + 1;
5190 diditb_offset = diditb_dim1 * diditb_dim2;
5191 diditb -= diditb_offset;
5192 soditb_dim1 = *nbpntu / 2;
5193 soditb_dim2 = *nbpntv / 2;
5194 soditb_offset = soditb_dim1 * (soditb_dim2 + 1) + 1;
5195 soditb -= soditb_offset;
5196 disotb_dim1 = *nbpntu / 2;
5197 disotb_dim2 = *nbpntv / 2;
5198 disotb_offset = disotb_dim1 * (disotb_dim2 + 1) + 1;
5199 disotb -= disotb_offset;
5200 sosotb_dim1 = *nbpntu / 2 + 1;
5201 sosotb_dim2 = *nbpntv / 2 + 1;
5202 sosotb_offset = sosotb_dim1 * sosotb_dim2;
5203 sosotb -= sosotb_offset;
5208 ldbg = AdvApp2Var_SysBase::mnfndeb_() >= 3;
5210 AdvApp2Var_SysBase::mgenmsg_("MMA2DS1", 7L);
5213 if (*isofav < 1 || *isofav > 2) {
5219 /* **********************************************************************
5221 /* --------- Discretization by U on the roots of the polynom of ------ */
5222 /* --------------- Legendre of degree NBPNTU, iso-V by iso-V --------- */
5223 /* **********************************************************************
5227 mma2ds2_(ndimen, &uintfn[1], &vintfn[1], foncnp, nbpntu, nbpntv, &
5228 urootb[1], &vrootb[1], &iuouv, &sosotb[sosotb_offset], &
5229 disotb[disotb_offset], &soditb[soditb_offset], &diditb[
5230 diditb_offset], &fpntab[fpntab_offset], &ttable[1], iercod);
5232 /* ******************************************************************
5234 /* --------- Discretization by V on the roots of the polynom of ------ */
5235 /* --------------- Legendre of degree NBPNTV, iso-V by iso-V --------- */
5236 /* ******************************************************************
5240 /* --> Inversion of indices of tables */
5242 for (nd = 1; nd <= i__1; ++nd) {
5243 isz1 = *nbpntu / 2 + 1;
5244 isz2 = *nbpntv / 2 + 1;
5245 AdvApp2Var_MathBase::mmfmtb1_(&isz1, &sosotb[nd * sosotb_dim2 * sosotb_dim1], &isz1, &
5246 isz2, &isz2, &sosotb[nd * sosotb_dim2 * sosotb_dim1], &
5247 ibid1, &ibid2, iercod);
5251 AdvApp2Var_MathBase::mmfmtb1_(&isz1, &diditb[nd * diditb_dim2 * diditb_dim1], &isz1, &
5252 isz2, &isz2, &diditb[nd * diditb_dim2 * diditb_dim1], &
5253 ibid1, &ibid2, iercod);
5259 AdvApp2Var_MathBase::mmfmtb1_(&isz1, &soditb[(nd * soditb_dim2 + 1) * soditb_dim1 + 1],
5260 &isz1, &isz2, &isz2, &soditb[(nd * soditb_dim2 + 1) *
5261 soditb_dim1 + 1], &ibid1, &ibid2, iercod);
5265 AdvApp2Var_MathBase::mmfmtb1_(&isz1, &disotb[(nd * disotb_dim2 + 1) * disotb_dim1 + 1],
5266 &isz1, &isz2, &isz2, &disotb[(nd * disotb_dim2 + 1) *
5267 disotb_dim1 + 1], &ibid1, &ibid2, iercod);
5274 mma2ds2_(ndimen, &vintfn[1], &uintfn[1], foncnp, nbpntv, nbpntu, &
5275 vrootb[1], &urootb[1], &iuouv, &sosotb[sosotb_offset], &
5276 soditb[soditb_offset], &disotb[disotb_offset], &diditb[
5277 diditb_offset], &fpntab[fpntab_offset], &ttable[1], iercod);
5278 /* --> Inversion of indices of tables */
5280 for (nd = 1; nd <= i__1; ++nd) {
5281 isz1 = *nbpntv / 2 + 1;
5282 isz2 = *nbpntu / 2 + 1;
5283 AdvApp2Var_MathBase::mmfmtb1_(&isz1, &sosotb[nd * sosotb_dim2 * sosotb_dim1], &isz1, &
5284 isz2, &isz2, &sosotb[nd * sosotb_dim2 * sosotb_dim1], &
5285 ibid1, &ibid2, iercod);
5289 AdvApp2Var_MathBase::mmfmtb1_(&isz1, &diditb[nd * diditb_dim2 * diditb_dim1], &isz1, &
5290 isz2, &isz2, &diditb[nd * diditb_dim2 * diditb_dim1], &
5291 ibid1, &ibid2, iercod);
5297 AdvApp2Var_MathBase::mmfmtb1_(&isz1, &soditb[(nd * soditb_dim2 + 1) * soditb_dim1 + 1],
5298 &isz1, &isz2, &isz2, &soditb[(nd * soditb_dim2 + 1) *
5299 soditb_dim1 + 1], &ibid1, &ibid2, iercod);
5303 AdvApp2Var_MathBase::mmfmtb1_(&isz1, &disotb[(nd * disotb_dim2 + 1) * disotb_dim1 + 1],
5304 &isz1, &isz2, &isz2, &disotb[(nd * disotb_dim2 + 1) *
5305 disotb_dim1 + 1], &ibid1, &ibid2, iercod);
5313 /* ------------------------------ The end -------------------------------
5319 AdvApp2Var_SysBase::maermsg_("MMA2DS1", iercod, 7L);
5322 AdvApp2Var_SysBase::mgsomsg_("MMA2DS1", 7L);
5327 //=======================================================================
5328 //function : mma2ds2_
5330 //=======================================================================
5331 int mma2ds2_(integer *ndimen,
5334 const AdvApp2Var_EvaluatorFunc2Var& foncnp,
5349 static integer c__0 = 0;
5350 /* System generated locals */
5351 integer sosotb_dim1, sosotb_dim2, sosotb_offset, disotb_dim1, disotb_dim2,
5352 disotb_offset, soditb_dim1, soditb_dim2, soditb_offset,
5353 diditb_dim1, diditb_dim2, diditb_offset, fpntab_dim1,
5354 fpntab_offset, i__1, i__2, i__3;
5356 /* Local variables */
5357 static integer jdec;
5358 static logical ldbg;
5359 static doublereal alinu, blinu, alinv, blinv, tcons;
5360 static doublereal dbfn1[2], dbfn2[2];
5361 static integer nuroo, nvroo, id, iu, iv;
5362 static doublereal um, up;
5365 /* **********************************************************************
5370 /* Discretization of function F(u,v) on the roots of polynoms of Legendre. */
5374 /* FONCTION&,DISCRETISATION,&POINT */
5376 /* INPUT ARGUMENTS : */
5377 /* ------------------ */
5378 /* NDIMEN: Dimension of the space. */
5379 /* UINTFN: Limits of the interval of definition by u of the function */
5380 /* to be processed: (UINTFN(1),UINTFN(2)). */
5381 /* VINTFN: Limits of the interval of definition by v of the function */
5382 /* to be processed: (VINTFN(1),VINTFN(2)). */
5383 /* FONCNP: The NAME of the non-polynomial function to be processed. */
5384 /* NBPNTU: The degree of Legendre polynom on the roots which of */
5385 /* FONCNP is discretized by u. */
5386 /* NBPNTV: The degree of Legendre polynom on the roots which of */
5387 /* FONCNP is discretized by v. */
5388 /* UROOTB: Table of STRICTLY POSITIVE roots of the polynom */
5389 /* of Legendre of degree NBPNTU defined on (-1,1). */
5390 /* VROOTB: Table of STRICTLY POSITIVE roots of the polynom */
5391 /* of Legendre of degree NBPNTV defined on (-1,1). */
5392 /* IIUOUV: Shows the type of iso of F(u,v) tom be extracted to improve the */
5393 /* rapidity of calculation (has no influence on the form of result) */
5394 /* = 1, shows that it is necessary to calculate the points of F(u,v) */
5395 /* with fixed u (so with NBPNTV values different from v). */
5396 /* = 2, shows that it is necessary to calculate the points of F(u,v) */
5397 /* with fixed v (so with NBPNTV values different from u). */
5398 /* SOSOTB: Preinitialized table (input/output argument). */
5399 /* DISOTB: Preinitialized table (input/output argument). */
5400 /* SODITB: Preinitialized table (input/output argument). */
5401 /* DIDITB: Preinitialized table (input/output argument). */
5403 /* OUTPUT ARGUMENTS : */
5404 /* ------------------- */
5405 /* SOSOTB: Table where the terms */
5406 /* F(ui,vj) + F(ui,-vj) + F(-ui,vj) + F(-ui,-vj) */
5407 /* are added with ui and vj positive roots of Legendre polynom */
5408 /* of degree NBPNTU and NBPNTV respectively. */
5409 /* DISOTB: Table where the terms */
5410 /* F(ui,vj) + F(ui,-vj) - F(-ui,vj) - F(-ui,-vj) */
5411 /* are added with ui and vj positive roots of Legendre polynom */
5412 /* of degree NBPNTU and NBPNTV respectively. */
5413 /* SODITB: Table where the terms */
5414 /* F(ui,vj) - F(ui,-vj) + F(-ui,vj) - F(-ui,-vj) */
5415 /* are added with ui and vj positive roots of Legendre polynom */
5416 /* of degree NBPNTU and NBPNTV respectively. */
5417 /* DIDITB: Table where the terms */
5418 /* F(ui,vj) - F(ui,-vj) - F(-ui,vj) + F(-ui,-vj) */
5419 /* are added with ui and vj positive roots of Legendre polynom */
5420 /* of degree NBPNTU and NBPNTV respectively. */
5421 /* FPNTAB: Auxiliary table. */
5422 /* TTABLE: Auxiliary table. */
5423 /* IERCOD: Error code >100 Pb in the evaluation of FONCNP, */
5424 /* the returned error code is equal to error code of FONCNP + 100. */
5426 /* COMMONS USED : */
5427 /* ---------------- */
5429 /* REFERENCES CALLED : */
5430 /* --------------------- */
5432 /* DESCRIPTION/NOTES/LIMITATIONS : */
5433 /* ----------------------------------- */
5434 /* --> The external function created by the caller of MA2F1K, MA2FDK */
5435 /* where MA2FXK should be in the following form : */
5436 /* SUBROUTINE FONCNP(NDIMEN,UINTFN,VINTFN,IIIUOUV,TCONST,NBPTAB */
5437 /* ,TTABLE,IDERIU,IDERIV,PPNTAB,IERCOD) */
5438 /* with the following input arguments : */
5439 /* - NDIMEN is integer defined as the sum of dimensions of */
5440 /* sub-spaces (i.e. total dimension of the problem). */
5441 /* - UINTFN(2) is a table of 2 reals containing the interval */
5442 /* by u where the function to be approximated is defined */
5443 /* (so it is equal to UIFONC). */
5444 /* - VINTFN(2) is a table of 2 reals containing the interval */
5445 /* by v where the function to be approximated is defined */
5446 /* (so it is equal to VIFONC). */
5447 /* - IIIUOUV, is 1 if it is necessary to calculate points with constant u, */
5448 /* is 2 if it is necessary to calculate points with constant v. */
5449 /* Any other value is an error. */
5450 /* - TCONST, real, value of the fixed parameter. Takes values */
5451 /* in (UIFONC(1),UIFONC(2)) if ISOFAV = 1 or */
5452 /* ins (VIFONC(1),VIFONC(2)) if ISOFAV = 2. */
5453 /* - NBPTAB, integer. Shows the number of points to be calculated. */
5454 /* - TTABLE, a table of reals NBPTAB. These are the values of */
5455 /* 'free' parameter of discretization (v if IIIUOUV=1, */
5456 /* u if IIIUOUV=2). */
5457 /* - IDERIU, integer, takes values between 0 (position) */
5458 /* and IORDRE(1) (partial derivative of the function by u */
5459 /* of order IORDRE(1) if IORDRE(1) > 0). */
5460 /* - IDERIV, integer, takes values between 0 (position) */
5461 /* and IORDRE(2) (partial derivative of the function by v */
5462 /* of order IORDRE(2) if IORDRE(2) > 0). */
5463 /* If IDERIU=i and IDERIV=j, FONCNP should calculate the */
5464 /* points of the derivative : */
5471 /* and the output arguments aret : */
5472 /* - FPNTAB(NDIMEN,NBPTAB) contains, at output, the table of */
5473 /* NBPTAB points calculated in FONCNP. */
5474 /* - IERCOD is, at output the error code of FONCNP. This code */
5475 /* (integer) should be strictly positive if there is a problem. */
5477 /* The input arguments SHOULD NOT be modified under FONCNP.
5480 /* -->As FONCNP is not forcedly defined in (-1,1)*(-1,1), the */
5481 /* values of UROOTB and VROOTB are consequently modified. */
5483 /* -->The results of discretisation are ranked in 4 tables */
5484 /* SOSOTB, DISOTB, SODITB and DIDITB to earn time */
5485 /* during the calculation of coefficients of the polynom of approximation. */
5487 /* When NBPNTU is uneven : */
5488 /* table SOSOTB(0,j) contains F(0,vj) + F(0,-vj), */
5489 /* table DIDITB(0,j) contains F(0,vj) - F(0,-vj), */
5490 /* When NBPNTV is uneven : */
5491 /* table SOSOTB(i,0) contains F(ui,0) + F(-ui,0), */
5492 /* table DIDITB(i,0) contains F(ui,0) - F(-ui,0), */
5493 /* When NBPNTU and NBPNTV are uneven : */
5494 /* term SOSOTB(0,0) contains F(0,0). */
5496 /* ATTENTION: These 4 tables are filled by varying the */
5497 /* 1st index first. So, the discretizations */
5498 /* of F(...,t) (for IIUOUV = 2) or of F(t,...) (IIUOUV = 1) */
5499 /* are stored in SOSOTB(...,t), SODITB(...,t), etc... */
5500 /* (this allows to gain important time). */
5501 /* It is required that the caller, in case of IIUOUV=1, */
5502 /* invert the roles of u and v, of SODITB and DISOTB BEFORE the */
5505 /* **********************************************************************
5508 /* Name of the routine */
5510 /* --> Indices of loops. */
5512 /* --------------------------- Initialization --------------------------
5515 /* Parameter adjustments */
5519 fpntab_dim1 = *ndimen;
5520 fpntab_offset = fpntab_dim1 + 1;
5521 fpntab -= fpntab_offset;
5523 diditb_dim1 = *nbpntu / 2 + 1;
5524 diditb_dim2 = *nbpntv / 2 + 1;
5525 diditb_offset = diditb_dim1 * diditb_dim2;
5526 diditb -= diditb_offset;
5527 soditb_dim1 = *nbpntu / 2;
5528 soditb_dim2 = *nbpntv / 2;
5529 soditb_offset = soditb_dim1 * (soditb_dim2 + 1) + 1;
5530 soditb -= soditb_offset;
5531 disotb_dim1 = *nbpntu / 2;
5532 disotb_dim2 = *nbpntv / 2;
5533 disotb_offset = disotb_dim1 * (disotb_dim2 + 1) + 1;
5534 disotb -= disotb_offset;
5535 sosotb_dim1 = *nbpntu / 2 + 1;
5536 sosotb_dim2 = *nbpntv / 2 + 1;
5537 sosotb_offset = sosotb_dim1 * sosotb_dim2;
5538 sosotb -= sosotb_offset;
5542 ldbg = AdvApp2Var_SysBase::mnfndeb_() >= 3;
5544 AdvApp2Var_SysBase::mgenmsg_("MMA2DS2", 7L);
5548 alinu = (uintfn[2] - uintfn[1]) / 2.;
5549 blinu = (uintfn[2] + uintfn[1]) / 2.;
5550 alinv = (vintfn[2] - vintfn[1]) / 2.;
5551 blinv = (vintfn[2] + vintfn[1]) / 2.;
5554 dbfn1[0] = vintfn[1];
5555 dbfn1[1] = vintfn[2];
5556 dbfn2[0] = uintfn[1];
5557 dbfn2[1] = uintfn[2];
5559 dbfn1[0] = uintfn[1];
5560 dbfn1[1] = uintfn[2];
5561 dbfn2[0] = vintfn[1];
5562 dbfn2[1] = vintfn[2];
5565 /* **********************************************************************
5567 /* -------- Discretization by U on the roots of Legendre polynom -------- */
5568 /* ---------------- of degree NBPNTU, with Vj fixed -------------------- */
5569 /* **********************************************************************
5572 nuroo = *nbpntu / 2;
5573 nvroo = *nbpntv / 2;
5574 jdec = (*nbpntu + 1) / 2;
5576 /* ----------- Loading of parameters of discretization by U ------------- */
5579 for (iu = 1; iu <= i__1; ++iu) {
5580 ttable[iu] = blinu + alinu * urootb[iu];
5584 /* -------------- For Vj fixed, negative root of Legendre ------------- */
5587 for (iv = 1; iv <= i__1; ++iv) {
5588 tcons = blinv + alinv * vrootb[iv];
5589 foncnp.Evaluate (ndimen, dbfn1, dbfn2, iiuouv, &tcons, nbpntu,
5590 &ttable[1], &c__0, &c__0, &fpntab[fpntab_offset], iercod);
5595 for (id = 1; id <= i__2; ++id) {
5597 for (iu = 1; iu <= i__3; ++iu) {
5598 up = fpntab[id + (iu + jdec) * fpntab_dim1];
5599 um = fpntab[id + (nuroo - iu + 1) * fpntab_dim1];
5600 sosotb[iu + (nvroo - iv + 1 + id * sosotb_dim2) * sosotb_dim1]
5601 = sosotb[iu + (nvroo - iv + 1 + id * sosotb_dim2) *
5602 sosotb_dim1] + up + um;
5603 disotb[iu + (nvroo - iv + 1 + id * disotb_dim2) * disotb_dim1]
5604 = disotb[iu + (nvroo - iv + 1 + id * disotb_dim2) *
5605 disotb_dim1] + up - um;
5606 soditb[iu + (nvroo - iv + 1 + id * soditb_dim2) * soditb_dim1]
5607 = soditb[iu + (nvroo - iv + 1 + id * soditb_dim2) *
5608 soditb_dim1] - up - um;
5609 diditb[iu + (nvroo - iv + 1 + id * diditb_dim2) * diditb_dim1]
5610 = diditb[iu + (nvroo - iv + 1 + id * diditb_dim2) *
5611 diditb_dim1] - up + um;
5614 if (*nbpntu % 2 != 0) {
5615 up = fpntab[id + jdec * fpntab_dim1];
5616 sosotb[(nvroo - iv + 1 + id * sosotb_dim2) * sosotb_dim1] +=
5618 diditb[(nvroo - iv + 1 + id * diditb_dim2) * diditb_dim1] -=
5626 /* --------- For Vj = 0 (uneven NBPNTV), discretization by U ----------- */
5628 if (*nbpntv % 2 != 0) {
5630 foncnp.Evaluate (ndimen, dbfn1, dbfn2, iiuouv, &tcons, nbpntu,
5631 &ttable[1], &c__0, &c__0, &fpntab[fpntab_offset], iercod);
5636 for (id = 1; id <= i__1; ++id) {
5638 for (iu = 1; iu <= i__2; ++iu) {
5639 up = fpntab[id + (jdec + iu) * fpntab_dim1];
5640 um = fpntab[id + (nuroo - iu + 1) * fpntab_dim1];
5641 sosotb[iu + id * sosotb_dim2 * sosotb_dim1] = sosotb[iu + id *
5642 sosotb_dim2 * sosotb_dim1] + up + um;
5643 diditb[iu + id * diditb_dim2 * diditb_dim1] = diditb[iu + id *
5644 diditb_dim2 * diditb_dim1] + up - um;
5647 if (*nbpntu % 2 != 0) {
5648 up = fpntab[id + jdec * fpntab_dim1];
5649 sosotb[id * sosotb_dim2 * sosotb_dim1] += up;
5655 /* -------------- For Vj fixed, positive root of Legendre ------------- */
5658 for (iv = 1; iv <= i__1; ++iv) {
5659 tcons = alinv * vrootb[(*nbpntv + 1) / 2 + iv] + blinv;
5660 foncnp.Evaluate (ndimen, dbfn1, dbfn2, iiuouv, &tcons, nbpntu,
5661 &ttable[1], &c__0, &c__0, &fpntab[fpntab_offset], iercod);
5666 for (id = 1; id <= i__2; ++id) {
5668 for (iu = 1; iu <= i__3; ++iu) {
5669 up = fpntab[id + (iu + jdec) * fpntab_dim1];
5670 um = fpntab[id + (nuroo - iu + 1) * fpntab_dim1];
5671 sosotb[iu + (iv + id * sosotb_dim2) * sosotb_dim1] = sosotb[
5672 iu + (iv + id * sosotb_dim2) * sosotb_dim1] + up + um;
5673 disotb[iu + (iv + id * disotb_dim2) * disotb_dim1] = disotb[
5674 iu + (iv + id * disotb_dim2) * disotb_dim1] + up - um;
5675 soditb[iu + (iv + id * soditb_dim2) * soditb_dim1] = soditb[
5676 iu + (iv + id * soditb_dim2) * soditb_dim1] + up + um;
5677 diditb[iu + (iv + id * diditb_dim2) * diditb_dim1] = diditb[
5678 iu + (iv + id * diditb_dim2) * diditb_dim1] + up - um;
5681 if (*nbpntu % 2 != 0) {
5682 up = fpntab[id + jdec * fpntab_dim1];
5683 sosotb[(iv + id * sosotb_dim2) * sosotb_dim1] += up;
5684 diditb[(iv + id * diditb_dim2) * diditb_dim1] += up;
5691 /* ------------------------------ The end -------------------------------
5697 AdvApp2Var_SysBase::maermsg_("MMA2DS2", iercod, 7L);
5700 AdvApp2Var_SysBase::mgsomsg_("MMA2DS2", 7L);
5705 //=======================================================================
5706 //function : mma2er1_
5708 //=======================================================================
5709 int mma2er1_(integer *ndjacu,
5725 /* System generated locals */
5726 integer patjac_dim1, patjac_dim2, patjac_offset, i__1, i__2, i__3;
5729 /* Local variables */
5730 static logical ldbg;
5731 static integer minu, minv;
5732 static doublereal vaux[2];
5733 static integer ii, nd, jj;
5734 static doublereal bid0, bid1;
5736 /* **********************************************************************
5741 /* Calculate max approximation error done when */
5742 /* the coefficients of PATJAC such that the degree by U varies between */
5743 /* MINDGU and MAXDGU and the degree by V varies between MINDGV and MAXDGV are removed. */
5747 /* TOUS,AB_SPECIFI:: CARREAU&,CALCUL,&ERREUR */
5749 /* INPUT ARGUMENTS : */
5750 /* ------------------ */
5751 /* NDJACU: Dimension by U of table PATJAC. */
5752 /* NDJACV: Dimension by V of table PATJAC. */
5753 /* NDIMEN: Dimension of the space. */
5754 /* MINDGU: Lower limit of index by U of coeff. of PATJAC to be taken into account. */
5755 /* MAXDGU: Upper limit of index by U of coeff. of PATJAC to be taken into account. */
5756 /* MINDGV: Lower limit of index by V of coeff. of PATJAC to be taken into account. */
5757 /* MAXDGV: Upper limit of index by V of coeff. of PATJAC to be taken into account. */
5758 /* IORDRU: Order of continuity by U provided by square PATJAC (from -1 to 2) */
5759 /* IORDRV: Order of continuity by U provided by square PATJAC (from -1 to 2) */
5760 /* XMAXJU: Maximum value of Jacobi polynoms of order IORDRU, */
5761 /* from degree 0 to MAXDGU - 2*(IORDU+1) */
5762 /* XMAXJV: Maximum value of Jacobi polynoms of order IORDRV, */
5763 /* from degree 0 to MAXDGV - 2*(IORDV+1) */
5764 /* PATJAC: Table of coeff. of square of approximation with */
5765 /* constraints of order IORDRU by U and IORDRV by V. */
5766 /* VECERR: Auxiliary vector. */
5767 /* ERREUR: MAX Error commited during removal of ALREADY CALCULATED coeff of PATJAC */
5769 /* OUTPUT ARGUMENTS : */
5770 /* ------------------- */
5771 /* ERREUR: MAX Error commited during removal of coeff of PATJAC */
5772 /* of indices from MINDGU to MAXDGU by U and from MINDGV to MAXDGV by V */
5773 /* THEN the already calculated error. */
5775 /* COMMONS USED : */
5776 /* ---------------- */
5778 /* REFERENCES CALLED : */
5779 /* --------------------- */
5781 /* DESCRIPTION/NOTES/LIMITATIONS : */
5782 /* ----------------------------------- */
5783 /* Table PATJAC is the place of storage of coeff. Cij of the square of */
5784 /* approximation of F(U,V). The indices i and j show the degree */
5785 /* by U and by V of base polynoms. These polynoms have the form: */
5787 /* ((1 - U*U)**(IORDRU+1)).J(i-2*(IORDRU+1)(U), where */
5789 /* polynom J(i-2*(IORDU+1)(U) is the Jacobi polynom of order */
5790 /* IORDRU+1 (the same by V by replacing U u V in the expression above). */
5792 /* The contribution to the error of term Cij when it is */
5793 /* removed from PATJAC is increased by: */
5795 /* DABS(Cij)*XMAXJU(i-2*(IORDRU+1))*XMAXJV(J-2*(IORDRV+1)) where */
5797 /* XMAXJU(i-2*(IORDRU+1) = ((1 - U*U)**(IORDRU+1)).J(i-2*(IORDRU+1)(U),
5799 /* XMAXJV(i-2*(IORDRV+1) = ((1 - V*V)**(IORDRV+1)).J(j-2*(IORDRV+1)(V).
5803 /* ***********************************************************************
5805 /* Name of the routine */
5808 /* ----------------------------- Initialisations ------------------------
5811 /* Parameter adjustments */
5813 patjac_dim1 = *ndjacu + 1;
5814 patjac_dim2 = *ndjacv + 1;
5815 patjac_offset = patjac_dim1 * patjac_dim2;
5816 patjac -= patjac_offset;
5819 ldbg = AdvApp2Var_SysBase::mnfndeb_() >= 3;
5821 AdvApp2Var_SysBase::mgenmsg_("MMA2ER1", 7L);
5824 minu = (*iordru + 1) << 1;
5825 minv = (*iordrv + 1) << 1;
5827 /* ------------------- Calculate the increment of the max error --------------- */
5828 /* ----- during the removal of the coeffs of indices from MINDGU to MAXDGU ---- */
5829 /* ---------------- by U and indices from MINDGV to MAXDGV by V --------------- */
5832 for (nd = 1; nd <= i__1; ++nd) {
5835 for (jj = *mindgv; jj <= i__2; ++jj) {
5838 for (ii = *mindgu; ii <= i__3; ++ii) {
5839 bid0 += (d__1 = patjac[ii + (jj + nd * patjac_dim2) *
5840 patjac_dim1], advapp_abs(d__1)) * xmaxju[ii - minu];
5843 bid1 = bid0 * xmaxjv[jj - minv] + bid1;
5851 /* ----------------------- Calculate the max error ----------------------*/
5853 bid1 = AdvApp2Var_MathBase::mzsnorm_(ndimen, &vecerr[1]);
5857 *erreur = AdvApp2Var_MathBase::mzsnorm_(&nd, vaux);
5859 /* ------------------------- The end ------------------------------------
5863 AdvApp2Var_SysBase::mgsomsg_("MMA2ER1", 7L);
5868 //=======================================================================
5869 //function : mma2er2_
5871 //=======================================================================
5872 int mma2er2_(integer *ndjacu,
5884 doublereal *epmscut,
5891 /* System generated locals */
5892 integer patjac_dim1, patjac_dim2, patjac_offset, i__1, i__2;
5895 /* Local variables */
5896 static logical ldbg;
5897 static doublereal vaux[2];
5898 static integer i2rdu, i2rdv;
5899 static doublereal errnu, errnv;
5900 static integer ii, nd, jj, nu, nv;
5901 static doublereal bid0, bid1;
5903 /* **********************************************************************
5908 /* Remove coefficients of PATJAC to obtain the minimum degree */
5909 /* by U and V checking the imposed tolerance. */
5913 /* TOUS,AB_SPECIFI:: CARREAU&,CALCUL,&ERREUR */
5915 /* INPUT ARGUMENTS : */
5916 /* ------------------ */
5917 /* NDJACU: Degree by U of table PATJAC. */
5918 /* NDJACV: Degree by V of table PATJAC. */
5919 /* NDIMEN: Dimension of the space. */
5920 /* MINDGU: Limit of index by U of coeff. of PATJAC to be PRESERVED (should be >=0). */
5921 /* MAXDGU: Upper limit of index by U of coeff. of PATJAC to be taken into account. */
5922 /* MINDGV: Limit of index by V of coeff. of PATJAC to be PRESERVED (should be >=0). */
5923 /* MAXDGV: Upper limit of index by V of coeff. of PATJAC to be taken into account. */
5924 /* IORDRU: Order of continuity by U provided by square PATJAC (from -1 to 2) */
5925 /* IORDRV: Order of continuity by U provided by square PATJAC (from -1 to 2) */
5926 /* XMAXJU: Maximum value of Jacobi polynoms of order IORDRU, */
5927 /* from degree 0 to MAXDGU - 2*(IORDU+1) */
5928 /* XMAXJV: Maximum value of Jacobi polynoms of order IORDRV, */
5929 /* from degree 0 to MAXDGV - 2*(IORDV+1) */
5930 /* PATJAC: Table of coeff. of square of approximation with */
5931 /* constraints of order IORDRU by U and IORDRV by V. */
5932 /* EPMSCUT: Tolerance of approximation. */
5933 /* VECERR: Auxiliary vector. */
5934 /* ERREUR: MAX Error commited ALREADY CALCULATED */
5936 /* OUTPUT ARGUMENTS : */
5937 /* ------------------- */
5938 /* ERREUR: MAX Error commited by preserving only coeff of PATJAC */
5939 /* of indices from 0 to NEWDGU by U and from 0 to NEWDGV by V */
5940 /* PLUS the already calculated error. */
5941 /* NEWDGU: Min. Degree by U such as the square of approximation */
5942 /* could check the tolerance. There is always NEWDGU >= MINDGU >= 0. */
5943 /* NEWDGV: Min. Degree by V such as the square of approximation */
5944 /* could check the tolerance. There is always NEWDGV >= MINDGV >= 0. */
5947 /* COMMONS USED : */
5948 /* ---------------- */
5950 /* REFERENCES CALLED : */
5951 /* --------------------- */
5953 /* DESCRIPTION/NOTES/LIMITATIONS : */
5954 /* ----------------------------------- */
5955 /* Table PATJAC is the place of storage of coeff. Cij of the square of */
5956 /* approximation of F(U,V). The indices i and j show the degree */
5957 /* by U and by V of base polynoms. These polynoms have the form: */
5959 /* ((1 - U*U)**(IORDRU+1)).J(i-2*(IORDRU+1)(U), where */
5961 /* polynom J(i-2*(IORDU+1)(U) is the Jacobi polynom of order */
5962 /* IORDRU+1 (the same by V by replacing U u V in the expression above). */
5964 /* The contribution to the error of term Cij when it is */
5965 /* removed from PATJAC is increased by: */
5967 /* DABS(Cij)*XMAXJU(i-2*(IORDRU+1))*XMAXJV(J-2*(IORDRV+1)) where */
5969 /* XMAXJU(i-2*(IORDRU+1) = ((1 - U*U)**(IORDRU+1)).J(i-2*(IORDRU+1)(U),
5971 /* XMAXJV(i-2*(IORDRV+1) = ((1 - V*V)**(IORDRV+1)).J(j-2*(IORDRV+1)(V).
5975 /* **********************************************************************
5977 /* Name of the routine */
5980 /* ----------------------------- Initialisations ------------------------
5983 /* Parameter adjustments */
5985 patjac_dim1 = *ndjacu + 1;
5986 patjac_dim2 = *ndjacv + 1;
5987 patjac_offset = patjac_dim1 * patjac_dim2;
5988 patjac -= patjac_offset;
5991 ldbg = AdvApp2Var_SysBase::mnfndeb_() >= 3;
5993 AdvApp2Var_SysBase::mgenmsg_("MMA2ER2", 7L);
5996 i2rdu = (*iordru + 1) << 1;
5997 i2rdv = (*iordrv + 1) << 1;
6001 /* **********************************************************************
6003 /* -------------------- Cutting of oefficients ------------------------
6005 /* **********************************************************************
6010 /* ------------------- Calculate the increment of max error --------------- */
6011 /* ----- during the removal of coeff. of indices from MINDGU to MAXDGU ------ */
6012 /* ---------------- by U, the degree by V is fixed to NV -----------------
6017 bid0 = xmaxjv[nv - i2rdv];
6019 for (nd = 1; nd <= i__1; ++nd) {
6022 for (ii = i2rdu; ii <= i__2; ++ii) {
6023 bid1 += (d__1 = patjac[ii + (nv + nd * patjac_dim2) *
6024 patjac_dim1], advapp_abs(d__1)) * xmaxju[ii - i2rdu] * bid0;
6031 vecerr[1] = *epmscut * 2;
6033 errnv = AdvApp2Var_MathBase::mzsnorm_(ndimen, &vecerr[1]);
6035 /* ------------------- Calculate the increment of max error --------------- */
6036 /* ----- during the removal of coeff. of indices from MINDGV to MAXDGV ------ */
6037 /* ---------------- by V, the degree by U is fixed to NU -----------------
6042 bid0 = xmaxju[nu - i2rdu];
6044 for (nd = 1; nd <= i__1; ++nd) {
6047 for (jj = i2rdv; jj <= i__2; ++jj) {
6048 bid1 += (d__1 = patjac[nu + (jj + nd * patjac_dim2) *
6049 patjac_dim1], advapp_abs(d__1)) * xmaxjv[jj - i2rdv] * bid0;
6056 vecerr[1] = *epmscut * 2;
6058 errnu = AdvApp2Var_MathBase::mzsnorm_(ndimen, &vecerr[1]);
6060 /* ----------------------- Calculate the max error ----------------------
6066 errnu = AdvApp2Var_MathBase::mzsnorm_(&nd, vaux);
6068 errnv = AdvApp2Var_MathBase::mzsnorm_(&nd, vaux);
6070 if (errnu > errnv) {
6071 if (errnv < *epmscut) {
6078 if (errnu < *epmscut) {
6088 /* -------------------------- Return the degrees -------------------
6092 *newdgu = advapp_max(nu,1);
6093 *newdgv = advapp_max(nv,1);
6095 /* ----------------------------------- The end --------------------------
6099 AdvApp2Var_SysBase::mgsomsg_("MMA2ER2", 7L);
6104 //=======================================================================
6105 //function : mma2fnc_
6107 //=======================================================================
6108 int AdvApp2Var_ApproxF2var::mma2fnc_(integer *ndimen,
6112 const AdvApp2Var_EvaluatorFunc2Var& foncnp,
6136 static integer c__8 = 8;
6138 /* System generated locals */
6139 integer courbe_dim1, courbe_dim2, courbe_offset, somtab_dim1, somtab_dim2,
6140 somtab_offset, diftab_dim1, diftab_dim2, diftab_offset,
6141 contr1_dim1, contr1_dim2, contr1_offset, contr2_dim1, contr2_dim2,
6142 contr2_offset, errmax_dim1, errmax_offset, errmoy_dim1,
6143 errmoy_offset, i__1;
6146 /* Local variables */
6147 static integer ideb;
6148 static doublereal tmil;
6149 static integer ideb1, ibid1, ibid2, ncfja, ndgre, ilong,
6151 static doublereal wrkar[1];
6152 static integer nupil;
6153 static long int iofwr;
6154 static doublereal uvpav[4] /* was [2][2] */;
6155 static integer nd, ii;
6158 static doublereal uv11[4] /* was [2][2] */;
6159 static integer ncb1;
6160 static doublereal eps3;
6161 static integer isz1, isz2, isz3, isz4, isz5;
6162 static long int ipt1, ipt2, ipt3, ipt4, ipt5,iptt, jptt;
6164 /* **********************************************************************
6169 /* Approximation of a limit of non polynomial function F(u,v) */
6170 /* (in the space of dimension NDIMEN) by SEVERAL */
6171 /* polynomial curves, by the method of least squares. The parameter of the function is preserved. */
6175 /* TOUS, AB_SPECIFI :: FONCTION&,EXTREMITE&, APPROXIMATION, &COURBE. */
6177 /* INPUT ARGUMENTS : */
6178 /* ----------------- */
6179 /* NDIMEN: Total Dimension of the space (sum of dimensions */
6180 /* of sub-spaces) */
6181 /* NBSESP: Number of "independent" sub-spaces. */
6182 /* NDIMSE: Table of dimensions of sub-spaces. */
6183 /* UVFONC: Limits of the interval (a,b)x(c,d) of definition of the */
6184 /* function to be approached by U (UVFONC(*,1) contains (a,b)) */
6185 /* and by V (UVFONC(*,2) contains (c,d)). */
6186 /* FONCNP: External function of position on the non polynomial function to be approached. */
6187 /* TCONST: Value of isoparameter of F(u,v) to be discretized. */
6188 /* ISOFAV: Type of chosen iso, = 1, shose that discretization is with u */
6189 /* fixed; = 2, shows that v is fixed. */
6190 /* NBROOT: Nb of points of discretisation of the iso, extremities not included. */
6191 /* ROOTLG: Table of roots of the polynom of Legendre defined on */
6192 /* (-1,1), of degree NBROOT. */
6193 /* IORDRE: Order of constraint at the extremities of the limit */
6194 /* -1 = no constraints, */
6195 /* 0 = constraints of passage to limits (i.e. C0), */
6196 /* 1 = C0 + constraints of 1st derivatives (i.e. C1), */
6197 /* 2 = C1 + constraints of 2nd derivatives (i.e. C2). */
6198 /* IDERIV: Order of derivative of the limit. */
6199 /* NDGJAC: Degree of serial development to be used for calculation in */
6200 /* the Jacobi base. */
6201 /* NBCRMX: Max Nb of curves to be created. */
6202 /* NCFLIM: Max Nb of coeff of the polynomial curve */
6203 /* of approximation (should be above or equal to */
6204 /* 2*IORDRE+2 and below or equal to 50). */
6205 /* EPSAPR: Table of required errors of approximation */
6206 /* sub-space by sub-space. */
6208 /* OUTPUT ARGUMENTS : */
6209 /* ------------------- */
6210 /* NCOEFF: Number of significative coeff of calculated curves. */
6211 /* COURBE: Table of coeff. of calculated polynomial curves. */
6212 /* Should be dimensioned in (NCFLIM,NDIMEN,NBCRMX). */
6213 /* These curves are ALWAYS parametrized in (-1,1). */
6214 /* NBCRBE: Nb of calculated curves. */
6215 /* SOMTAB: For F defined on (-1,1) (otherwise rescale the */
6216 /* parameters), this is the table of sums F(u,vj) + F(u,-vj)
6218 /* if ISOFAV = 1 (and IDERIV=0, otherwise the derivatives */
6219 /* by u of order IDERIV are taken) or sumes F(ui,v) + F(-ui,v) if */
6220 /* ISOFAV = 2 (and IDERIV=0, otherwise the derivatives by */
6221 /* v of order IDERIV are taken). */
6222 /* DIFTAB: For F defined on (-1,1) (otherwise rescale the */
6223 /* parameters), this is the table of sums F(u,vj) - F(u,-vj)
6225 /* if ISOFAV = 1 (and IDERIV=0, otherwise the derivatives */
6226 /* by u of order IDERIV are taken) or sumes F(ui,v) + F(-ui,v) if */
6227 /* ISOFAV = 2 (and IDERIV=0, otherwise the derivatives by */
6228 /* v of order IDERIV are taken). */
6229 /* CONTR1: Contains the coordinates of the left extremity of the iso */
6230 /* and of its derivatives till order IORDRE */
6231 /* CONTR2: Contains the coordinates of the right extremity of the iso */
6232 /* and of its derivatives till order IORDRE */
6233 /* TABDEC: Table of NBCRBE+1 parameters of cut of UVFONC(1:2,1)
6235 /* if ISOFAV=2, or of UVFONC(1:2,2) if ISOFAV=1. */
6236 /* ERRMAX: Table of MAX errors (sub-space by sub-space) */
6237 /* committed in the approximation of FONCNP by NBCRBE curves. */
6238 /* ERRMOY: Table of AVERAGE errors (sub-space by sub-space) */
6239 /* committed in the approximation of FONCNP by NBCRBE curves.
6240 /* IERCOD: Error code: */
6241 /* -1 = ERRMAX > EPSAPR for at least one sub-space. */
6242 /* (the resulting curves of at least mathematic degree NCFLIM-1 */
6243 /* are calculated). */
6244 /* 0 = Everything is ok. */
6245 /* 1 = Pb of incoherence of inputs. */
6246 /* 10 = Pb of calculation of the interpolation of constraints. */
6247 /* 13 = Pb in the dynamic allocation. */
6248 /* 33 = Pb in the data recuperation from block data */
6249 /* of coeff. of integration by GAUSS method. */
6250 /* >100 Pb in the evaluation of FONCNP, the returned error code */
6251 /* is equal to the error code of FONCNP + 100. */
6253 /* COMMONS USED : */
6254 /* ---------------- */
6256 /* REFERENCES CALLED : */
6257 /* ----------------------- */
6259 /* DESCRIPTION/NOTES/LIMITATIONS : */
6260 /* ----------------------------------- */
6261 /* --> The approximation part is done in the space of dimension */
6262 /* NDIMEN (the sum of dimensions of sub-spaces). For example : */
6263 /* If NBSESP=2 and NDIMSE(1)=3, NDIMSE(2)=2, there is smoothing with */
6264 /* NDIMEN=5. The result (in COURBE(NDIMEN,NCOEFF,i) ), will be */
6265 /* composed of the result of smoothing of 3D function in */
6266 /* COURBE(1:3,1:NCOEFF,i) and of smoothing of 2D function in */
6267 /* COURBE(4:5,1:NCOEFF,i). */
6269 /* --> Routine FONCNP should be declared EXTERNAL in the program */
6270 /* calling MMA2FNC. */
6272 /* --> Function FONCNP, declared externally, should be declared */
6273 /* IMPERATIVELY in form : */
6274 /* SUBROUTINE FONCNP(NDIMEN,UINTFN,VINTFN,IIUOUV,TCONST,NBPTAB */
6275 /* ,TTABLE,IDERIU,IDERIV,IERCOD) */
6276 /* where the input arguments are : */
6277 /* - NDIMEN is integer defined as the sum of dimensions of */
6278 /* sub-spaces (i.e. total dimension of the problem). */
6279 /* - UINTFN(2) is a table of 2 reals containing the interval */
6280 /* by u where the function to be approximated is defined */
6281 /* (so it is equal to UIFONC). */
6282 /* - VINTFN(2) is a table of 2 reals containing the interval */
6283 /* by v where the function to be approximated is defined */
6284 /* (so it is equal to VIFONC). */
6285 /* - IIUOUV, shows that the points to be calculated have a constant U */
6286 /* (IIUOUV=1) or a constant V (IIUOUV=2). */
6287 /* - TCONST, real, value of the fixed discretisation parameter. Takes values */
6288 /* in (UINTFN(1),UINTFN(2)) if IIUOUV=1, */
6289 /* or in (VINTFN(1),VINTFN(2)) if IIUOUV=2. */
6290 /* - NBPTAB, the nb of point of discretisation following the free variable */
6291 /* : V if IIUOUV=1 or U if IIUOUV = 2. */
6292 /* - TTABLE, Table of NBPTAB parametres of discretisation. . */
6293 /* - IDERIU, integer, takes values between 0 (position) */
6294 /* and IORDREU (partial derivative of the function by u */
6295 /* of order IORDREU if IORDREU > 0). */
6296 /* - IDERIV, integer, takes values between 0 (position) */
6297 /* and IORDREV (partial derivative of the function by v */
6298 /* of order IORDREV if IORDREV > 0). */
6299 /* and the output arguments are : */
6300 /* - FPNTAB(NDIMEN,NBPTAB) contains, at output, the table of */
6301 /* NBPTAB points calculated in FONCNP. */
6302 /* - IERCOD is, at output the error code of FONCNP. This code */
6303 /* (integer) should be strictly positive if there is a problem. */
6305 /* The input arguments SHOULD NOT BE modified under FONCNP.
6308 /* --> If IERCOD=-1, the required precision can't be reached (ERRMAX */
6309 /* is above EPSAPR on at least one sub-space), but
6311 /* one gives the best possible result for NCFLIM and EPSAPR */
6312 /* chosen by the user. In this case (and for IERCOD=0), there is a solution. */
6315 /* **********************************************************************
6317 /* Name of the routine */
6319 /* Parameter adjustments */
6324 errmoy_dim1 = *nbsesp;
6325 errmoy_offset = errmoy_dim1 + 1;
6326 errmoy -= errmoy_offset;
6327 errmax_dim1 = *nbsesp;
6328 errmax_offset = errmax_dim1 + 1;
6329 errmax -= errmax_offset;
6330 contr2_dim1 = *ndimen;
6331 contr2_dim2 = *iordre + 2;
6332 contr2_offset = contr2_dim1 * (contr2_dim2 + 1) + 1;
6333 contr2 -= contr2_offset;
6334 contr1_dim1 = *ndimen;
6335 contr1_dim2 = *iordre + 2;
6336 contr1_offset = contr1_dim1 * (contr1_dim2 + 1) + 1;
6337 contr1 -= contr1_offset;
6338 diftab_dim1 = *nbroot / 2 + 1;
6339 diftab_dim2 = *ndimen;
6340 diftab_offset = diftab_dim1 * (diftab_dim2 + 1);
6341 diftab -= diftab_offset;
6342 somtab_dim1 = *nbroot / 2 + 1;
6343 somtab_dim2 = *ndimen;
6344 somtab_offset = somtab_dim1 * (somtab_dim2 + 1);
6345 somtab -= somtab_offset;
6347 courbe_dim1 = *ncflim;
6348 courbe_dim2 = *ndimen;
6349 courbe_offset = courbe_dim1 * (courbe_dim2 + 1) + 1;
6350 courbe -= courbe_offset;
6353 ibb = AdvApp2Var_SysBase::mnfndeb_();
6355 AdvApp2Var_SysBase::mgenmsg_("MMA2FNC", 7L);
6360 /* ---------------- Set to zero the coefficients of CURVE --------------
6363 ilong = *ndimen * *ncflim * *nbcrmx;
6364 AdvApp2Var_SysBase::mvriraz_(&ilong, (char *)&courbe[courbe_offset]);
6366 /* **********************************************************************
6368 /* -------------------------- Checking of entries ------------------
6370 /* **********************************************************************
6373 AdvApp2Var_MathBase::mmveps3_(&eps3);
6374 if ((d__1 = uvfonc[4] - uvfonc[3], advapp_abs(d__1)) < eps3) {
6377 if ((d__1 = uvfonc[6] - uvfonc[5], advapp_abs(d__1)) < eps3) {
6386 /* ********************************************************************** */
6387 /* ------------- Preparation of parameters of discretisation ----------- */
6388 /* **********************************************************************
6391 /* -- Allocation of a table of parameters and points of discretisation -- */
6392 /* --> For the parameters of discretisation. */
6394 /* --> For the points of discretisation in MMA1FDI and MMA1CDI and
6396 /* the auxiliary curve for MMAPCMP */
6397 ibid1 = *ndimen * (*nbroot + 2);
6398 ibid2 = ((*iordre + 1) << 1) * *nbroot;
6399 isz2 = advapp_max(ibid1,ibid2);
6400 ibid1 = (((*ncflim - 1) / 2 + 1) << 1) * *ndimen;
6401 isz2 = advapp_max(ibid1,isz2);
6402 /* --> To return the polynoms of hermit. */
6403 isz3 = ((*iordre + 1) << 2) * (*iordre + 1);
6404 /* --> For the Gauss coeff. of integration. */
6405 isz4 = (*nbroot / 2 + 1) * (*ndgjac + 1 - ((*iordre + 1) << 1));
6406 /* --> For the coeff of the curve in the base of Jacobi */
6407 isz5 = (*ndgjac + 1) * *ndimen;
6409 ndwrk = isz1 + isz2 + isz3 + isz4 + isz5;
6410 AdvApp2Var_SysBase::mcrrqst_(&c__8, &ndwrk, wrkar, &iofwr, &ier);
6413 /* --> For the parameters of discretisation (NBROOT+2 extremities). */
6415 /* --> For the points of discretisation FPNTAB(NDIMEN,NBROOT+2), */
6416 /* FPNTAB(NBROOT,2*(IORDRE+1)) and for WRKAR of MMAPCMP. */
6418 /* --> For the polynoms of Hermit */
6420 /* --> For the Gauss coeff of integration. */
6422 /* --> For the curve in Jacobi. */
6425 /* ------------------ Initialisation of management of cuts ---------
6429 uvpav[0] = uvfonc[3];
6430 uvpav[1] = uvfonc[4];
6431 tabdec[0] = uvfonc[5];
6432 tabdec[1] = uvfonc[6];
6433 } else if (*isofav == 2) {
6434 tabdec[0] = uvfonc[3];
6435 tabdec[1] = uvfonc[4];
6436 uvpav[2] = uvfonc[5];
6437 uvpav[3] = uvfonc[6];
6445 /* **********************************************************************
6447 /* APPROXIMATION WITH CUTS */
6448 /* **********************************************************************
6452 /* --> When the top is reached, this is the end ! */
6453 if (nupil - *nbcrbe == 0) {
6458 uvpav[2] = tabdec[*nbcrbe];
6459 uvpav[3] = tabdec[*nbcrbe + 1];
6460 } else if (*isofav == 2) {
6461 uvpav[0] = tabdec[*nbcrbe];
6462 uvpav[1] = tabdec[*nbcrbe + 1];
6467 /* -------------------- Normalization of parameters -------------------- */
6469 mma1nop_(nbroot, &rootlg[1], uvpav, isofav, &wrkar[ipt1], &ier);
6474 /* -------------------- Discretisation of FONCNP ------------------------ */
6476 mma1fdi_(ndimen, uvpav, foncnp, isofav, tconst, nbroot, &wrkar[ipt1],
6477 iordre, ideriv, &wrkar[ipt2], &somtab[(ncb1 * somtab_dim2 + 1) *
6478 somtab_dim1], &diftab[(ncb1 * diftab_dim2 + 1) * diftab_dim1], &
6479 contr1[(ncb1 * contr1_dim2 + 1) * contr1_dim1 + 1], &contr2[(ncb1
6480 * contr2_dim2 + 1) * contr2_dim1 + 1], iercod);
6485 /* -----------Cut the discretisation of constraints ------------*/
6488 mma1cdi_(ndimen, nbroot, &rootlg[1], iordre, &contr1[(ncb1 *
6489 contr1_dim2 + 1) * contr1_dim1 + 1], &contr2[(ncb1 *
6490 contr2_dim2 + 1) * contr2_dim1 + 1], &somtab[(ncb1 *
6491 somtab_dim2 + 1) * somtab_dim1], &diftab[(ncb1 * diftab_dim2
6492 + 1) * diftab_dim1], &wrkar[ipt2], &wrkar[ipt3], &ier);
6498 /* **********************************************************************
6500 /* -------------------- Calculate the curve of approximation -------------
6502 /* **********************************************************************
6505 mma1jak_(ndimen, nbroot, iordre, ndgjac, &somtab[(ncb1 * somtab_dim2 + 1)
6506 * somtab_dim1], &diftab[(ncb1 * diftab_dim2 + 1) * diftab_dim1], &
6507 wrkar[ipt4], &wrkar[ipt5], &ier);
6512 /* **********************************************************************
6514 /* ---------------- Add polynom of interpolation -------------------
6516 /* **********************************************************************
6520 mma1cnt_(ndimen, iordre, &contr1[(ncb1 * contr1_dim2 + 1) *
6521 contr1_dim1 + 1], &contr2[(ncb1 * contr2_dim2 + 1) *
6522 contr2_dim1 + 1], &wrkar[ipt3], ndgjac, &wrkar[ipt5]);
6525 /* **********************************************************************
6527 /* --------------- Calculate Max and Average error ----------------------
6529 /* **********************************************************************
6532 mma1fer_(ndimen, nbsesp, &ndimse[1], iordre, ndgjac, &wrkar[ipt5], ncflim,
6533 &epsapr[1], &wrkar[ipt2], &errmax[ncb1 * errmax_dim1 + 1], &
6534 errmoy[ncb1 * errmoy_dim1 + 1], &ncoeff[ncb1], &ier);
6539 if (ier == 0 || (ier == -1 && nupil == *nbcrmx)) {
6541 /* ******************************************************************
6543 /* ----------------------- Compression du resultat ------------------
6545 /* ******************************************************************
6551 ncfja = *ndgjac + 1;
6552 /* -> Compression of result in WRKAR(IPT2) */
6555 AdvApp2Var_MathBase::mmapcmp_(ndimen,
6556 &ncfja, &ncoeff[ncb1], &wrkar[ipt5], &wrkar[ipt2]);
6558 AdvApp2Var_MathBase::mmapcmp_((integer*)ndimen,
6564 ilong = *ndimen * *ncflim;
6565 AdvApp2Var_SysBase::mvriraz_(&ilong, (char*)&wrkar[ipt5]);
6566 /* -> Passage to canonic base (-1,1) (result in WRKAR(IPT5)).
6568 ndgre = ncoeff[ncb1] - 1;
6570 for (nd = 1; nd <= i__1; ++nd) {
6571 iptt = ipt2 + ((nd - 1) << 1) * (ndgre / 2 + 1);
6572 jptt = ipt5 + (nd - 1) * ncoeff[ncb1];
6573 AdvApp2Var_MathBase::mmjacan_(iordre, &ndgre, &wrkar[iptt], &wrkar[jptt]);
6577 /* -> Store the calculated curve */
6579 AdvApp2Var_MathBase::mmfmca8_(&ncoeff[ncb1], ndimen, &ibid1, ncflim, ndimen, &ibid1, &
6580 wrkar[ipt5], &courbe[(ncb1 * courbe_dim2 + 1) * courbe_dim1 +
6583 /* -> Before normalization of constraints on (-1,1), recalculate */
6584 /* the true constraints. */
6586 for (ii = 0; ii <= i__1; ++ii) {
6587 mma1noc_(uv11, ndimen, &ii, &contr1[(ii + 1 + ncb1 * contr1_dim2)
6588 * contr1_dim1 + 1], uvpav, isofav, ideriv, &contr1[(ii +
6589 1 + ncb1 * contr1_dim2) * contr1_dim1 + 1]);
6590 mma1noc_(uv11, ndimen, &ii, &contr2[(ii + 1 + ncb1 * contr2_dim2)
6591 * contr2_dim1 + 1], uvpav, isofav, ideriv, &contr2[(ii +
6592 1 + ncb1 * contr2_dim2) * contr2_dim1 + 1]);
6596 ibid1 = (*nbroot / 2 + 1) * *ndimen;
6597 mma1noc_(uv11, &ibid1, &ii, &somtab[(ncb1 * somtab_dim2 + 1) *
6598 somtab_dim1], uvpav, isofav, ideriv, &somtab[(ncb1 *
6599 somtab_dim2 + 1) * somtab_dim1]);
6600 mma1noc_(uv11, &ibid1, &ii, &diftab[(ncb1 * diftab_dim2 + 1) *
6601 diftab_dim1], uvpav, isofav, ideriv, &diftab[(ncb1 *
6602 diftab_dim2 + 1) * diftab_dim1]);
6605 for (nd = 1; nd <= i__1; ++nd) {
6606 mma1noc_(uv11, &ncoeff[ncb1], &ii, &courbe[(nd + ncb1 *
6607 courbe_dim2) * courbe_dim1 + 1], uvpav, isofav, ideriv, &
6608 courbe[(nd + ncb1 * courbe_dim2) * courbe_dim1 + 1]);
6612 /* -> Update the nb of already created curves */
6615 /* -> ...otherwise try to cut the current interval in 2... */
6617 tmil = (tabdec[*nbcrbe + 1] + tabdec[*nbcrbe]) / 2.;
6620 ilong = (nupil - *nbcrbe) << 3;
6621 AdvApp2Var_SysBase::mcrfill_(&ilong, (char *)&tabdec[ideb],(char *)&tabdec[ideb1]);
6622 tabdec[ideb] = tmil;
6626 /* ---------- Make approximation of the rest -----------
6631 /* --------------------- Return code of error -----------------
6633 /* --> Pb with dynamic allocation */
6637 /* --> Inputs incoherent. */
6642 /* -------------------------- Dynamic desallocation -------------------
6647 AdvApp2Var_SysBase::mcrdelt_(&c__8, &ndwrk, wrkar, &iofwr, &ier);
6654 /* ------------------------------ The end -------------------------------
6659 AdvApp2Var_SysBase::maermsg_("MMA2FNC", iercod, 7L);
6662 AdvApp2Var_SysBase::mgsomsg_("MMA2FNC", 7L);
6667 //=======================================================================
6668 //function : mma2fx6_
6670 //=======================================================================
6671 int AdvApp2Var_ApproxF2var::mma2fx6_(integer *ncfmxu,
6688 /* System generated locals */
6689 integer epsfro_dim1, epsfro_offset, patcan_dim1, patcan_dim2, patcan_dim3,
6690 patcan_dim4, patcan_offset, errmax_dim1, errmax_dim2,
6691 errmax_offset, ncoefu_dim1, ncoefu_offset, ncoefv_dim1,
6692 ncoefv_offset, i__1, i__2, i__3, i__4, i__5;
6693 doublereal d__1, d__2;
6695 /* Local variables */
6696 static integer idim, ncfu, ncfv, id, ii, nd, jj, ku, kv, ns, ibb;
6697 static doublereal bid;
6698 static doublereal tol;
6700 /* **********************************************************************
6705 /* Reduction of degree when the squares are the squares of constraints. */
6709 /* TOUS,AB_SPECIFI::CARREAU&,REDUCTION,&CARREAU */
6711 /* INPUT ARGUMENTS : */
6712 /* ------------------ */
6713 /* NCFMXU: Max Nb of coeff by u of solution P(u,v) (table */
6714 /* PATCAN). This argument serves only to declare the size of this table. */
6715 /* NCFMXV: Max Nb of coeff by v of solution P(u,v) (table */
6716 /* PATCAN). This argument serves only to declare the size of this table. */
6717 /* NDIMEN: Total dimension of the space where the processed function */
6718 /* takes its values.(sum of dimensions of sub-spaces) */
6719 /* NBSESP: Nb of independent sub-spaces where the errors are measured. */
6720 /* NDIMSE: Table of dimensions of NBSESP sub-spaces. */
6721 /* NBUPAT: Nb of square solution by u. */
6722 /* NBVPAT: Nb of square solution by v. */
6723 /* IORDRU: Order of constraint imposed at the extremities of iso-V */
6724 /* = 0, the extremities of iso-V are calculated */
6725 /* = 1, additionally the 1st derivative in the direction of iso-V is calculated */
6726 /* = 2, additionally the 2nd derivative in the direction of iso-V is calculated */
6727 /* IORDRV: Ordre de contrainte impose aux extremites de l'iso-U */
6728 /* = 0, on calcule les extremites de l'iso-U. */
6729 /* = 1, additionally the 1st derivative in the direction of iso-U is calculated */
6730 /* = 2, additionally the 2nd derivative in the direction of iso-U is calculated */
6731 /* EPSAPR: Table of imposed precisions, sub-space by sub-space. */
6732 /* EPSFRO: Table of imposed precisions, sub-space by sub-space on the limits of squares. */
6733 /* PATCAN: Table of coeff. in the canonic base of squares P(u,v) calculated for (u,v) in (-1,1). */
6734 /* ERRMAX: Table of MAX errors (sub-space by sub-space) */
6735 /* committed in the approximation of F(u,v) by P(u,v). */
6736 /* NCOEFU: Table of Nb of significative coeffs. by u of calculated squares. */
6737 /* NCOEFV: Table of Nb of significative coeffs. by v of calculated squares. */
6739 /* OUTPUT ARGUMENTS : */
6740 /* ------------------- */
6741 /* NCOEFU: Table of Nb of significative coeffs. by u of calculated squares. */
6742 /* NCOEFV: Table of Nb of significative coeffs. by v of calculated squares. */
6744 /* COMMONS USED : */
6745 /* ---------------- */
6747 /* REFERENCES CALLED : */
6748 /* --------------------- */
6750 /* DESCRIPTION/NOTES/LIMITATIONS : */
6751 /* ------------------------------- */
6753 /* **********************************************************************
6756 /* Name of the routine */
6759 /* Parameter adjustments */
6760 epsfro_dim1 = *nbsesp;
6761 epsfro_offset = epsfro_dim1 * 5 + 1;
6762 epsfro -= epsfro_offset;
6765 ncoefv_dim1 = *nbupat;
6766 ncoefv_offset = ncoefv_dim1 + 1;
6767 ncoefv -= ncoefv_offset;
6768 ncoefu_dim1 = *nbupat;
6769 ncoefu_offset = ncoefu_dim1 + 1;
6770 ncoefu -= ncoefu_offset;
6771 errmax_dim1 = *nbsesp;
6772 errmax_dim2 = *nbupat;
6773 errmax_offset = errmax_dim1 * (errmax_dim2 + 1) + 1;
6774 errmax -= errmax_offset;
6775 patcan_dim1 = *ncfmxu;
6776 patcan_dim2 = *ncfmxv;
6777 patcan_dim3 = *ndimen;
6778 patcan_dim4 = *nbupat;
6779 patcan_offset = patcan_dim1 * (patcan_dim2 * (patcan_dim3 * (patcan_dim4
6781 patcan -= patcan_offset;
6784 ibb = AdvApp2Var_SysBase::mnfndeb_();
6786 AdvApp2Var_SysBase::mgenmsg_("MMA2FX6", 7L);
6791 for (jj = 1; jj <= i__1; ++jj) {
6793 for (ii = 1; ii <= i__2; ++ii) {
6794 ncfu = ncoefu[ii + jj * ncoefu_dim1];
6795 ncfv = ncoefv[ii + jj * ncoefv_dim1];
6797 /* ********************************************************************** */
6798 /* -------------------- Reduction of degree by U ------------------------- */
6799 /* ********************************************************************** */
6802 if (ncfu <= (*iordru + 1) << 1 && ncfu > 2) {
6806 for (ns = 1; ns <= i__3; ++ns) {
6809 d__1 = tol, d__2 = epsfro[ns + epsfro_dim1 * 9];
6810 tol = advapp_min(d__1,d__2);
6812 d__1 = tol, d__2 = epsfro[ns + epsfro_dim1 * 10];
6813 tol = advapp_min(d__1,d__2);
6815 d__1 = tol, d__2 = epsfro[ns + epsfro_dim1 * 11];
6816 tol = advapp_min(d__1,d__2);
6818 d__1 = tol, d__2 = epsfro[ns + epsfro_dim1 * 12];
6819 tol = advapp_min(d__1,d__2);
6820 if (ii == 1 || ii == *nbupat || jj == 1 || jj == *nbvpat)
6823 d__1 = tol, d__2 = epsfro[ns + epsfro_dim1 * 5];
6824 tol = advapp_min(d__1,d__2);
6826 d__1 = tol, d__2 = epsfro[ns + epsfro_dim1 * 6];
6827 tol = advapp_min(d__1,d__2);
6829 d__1 = tol, d__2 = epsfro[ns + epsfro_dim1 * 7];
6830 tol = advapp_min(d__1,d__2);
6832 d__1 = tol, d__2 = epsfro[ns + (epsfro_dim1 << 3)];
6833 tol = advapp_min(d__1,d__2);
6838 for (nd = 1; nd <= i__4; ++nd) {
6841 for (kv = 1; kv <= i__5; ++kv) {
6842 bid += (d__1 = patcan[ncfu + (kv + (id + (ii + jj
6843 * patcan_dim4) * patcan_dim3) *
6844 patcan_dim2) * patcan_dim1], advapp_abs(d__1));
6850 if (bid > tol * 1e-6 || bid > errmax[ns + (ii + jj *
6851 errmax_dim2) * errmax_dim1]) {
6862 /* ********************************************************************** */
6863 /* -------------------- Reduction of degree by V ------------------------- */
6864 /* ********************************************************************** */
6867 if (ncfv <= (*iordrv + 1) << 1 && ncfv > 2) {
6871 for (ns = 1; ns <= i__3; ++ns) {
6874 d__1 = tol, d__2 = epsfro[ns + epsfro_dim1 * 9];
6875 tol = advapp_min(d__1,d__2);
6877 d__1 = tol, d__2 = epsfro[ns + epsfro_dim1 * 10];
6878 tol = advapp_min(d__1,d__2);
6880 d__1 = tol, d__2 = epsfro[ns + epsfro_dim1 * 11];
6881 tol = advapp_min(d__1,d__2);
6883 d__1 = tol, d__2 = epsfro[ns + epsfro_dim1 * 12];
6884 tol = advapp_min(d__1,d__2);
6885 if (ii == 1 || ii == *nbupat || jj == 1 || jj == *nbvpat)
6888 d__1 = tol, d__2 = epsfro[ns + epsfro_dim1 * 5];
6889 tol = advapp_min(d__1,d__2);
6891 d__1 = tol, d__2 = epsfro[ns + epsfro_dim1 * 6];
6892 tol = advapp_min(d__1,d__2);
6894 d__1 = tol, d__2 = epsfro[ns + epsfro_dim1 * 7];
6895 tol = advapp_min(d__1,d__2);
6897 d__1 = tol, d__2 = epsfro[ns + (epsfro_dim1 << 3)];
6898 tol = advapp_min(d__1,d__2);
6903 for (nd = 1; nd <= i__4; ++nd) {
6906 for (ku = 1; ku <= i__5; ++ku) {
6907 bid += (d__1 = patcan[ku + (ncfv + (id + (ii + jj
6908 * patcan_dim4) * patcan_dim3) *
6909 patcan_dim2) * patcan_dim1], advapp_abs(d__1));
6915 if (bid > tol * 1e-6 || bid > errmax[ns + (ii + jj *
6916 errmax_dim2) * errmax_dim1]) {
6927 /* --- Return the nbs of coeff. and pass to the next square --- */
6930 ncoefu[ii + jj * ncoefu_dim1] = advapp_max(ncfu,2);
6931 ncoefv[ii + jj * ncoefv_dim1] = advapp_max(ncfv,2);
6937 /* ------------------------------ The End -------------------------------
6941 AdvApp2Var_SysBase::mgsomsg_("MMA2FX6", 7L);
6947 //=======================================================================
6948 //function : mma2jmx_
6950 //=======================================================================
6951 int AdvApp2Var_ApproxF2var::mma2jmx_(integer *ndgjac,
6955 /* Initialized data */
6957 static doublereal xmax2[57] = { .9682458365518542212948163499456,
6958 .986013297183269340427888048593603,
6959 1.07810420343739860362585159028115,
6960 1.17325804490920057010925920756025,
6961 1.26476561266905634732910520370741,
6962 1.35169950227289626684434056681946,
6963 1.43424378958284137759129885012494,
6964 1.51281316274895465689402798226634,
6965 1.5878364329591908800533936587012,
6966 1.65970112228228167018443636171226,
6967 1.72874345388622461848433443013543,
6968 1.7952515611463877544077632304216,
6969 1.85947199025328260370244491818047,
6970 1.92161634324190018916351663207101,
6971 1.98186713586472025397859895825157,
6972 2.04038269834980146276967984252188,
6973 2.09730119173852573441223706382076,
6974 2.15274387655763462685970799663412,
6975 2.20681777186342079455059961912859,
6976 2.25961782459354604684402726624239,
6977 2.31122868752403808176824020121524,
6978 2.36172618435386566570998793688131,
6979 2.41117852396114589446497298177554,
6980 2.45964731268663657873849811095449,
6981 2.50718840313973523778244737914028,
6982 2.55385260994795361951813645784034,
6983 2.59968631659221867834697883938297,
6984 2.64473199258285846332860663371298,
6985 2.68902863641518586789566216064557,
6986 2.73261215675199397407027673053895,
6987 2.77551570192374483822124304745691,
6988 2.8177699459714315371037628127545,
6989 2.85940333797200948896046563785957,
6990 2.90044232019793636101516293333324,
6991 2.94091151970640874812265419871976,
6992 2.98083391718088702956696303389061,
6993 3.02023099621926980436221568258656,
6994 3.05912287574998661724731962377847,
6995 3.09752842783622025614245706196447,
6996 3.13546538278134559341444834866301,
6997 3.17295042316122606504398054547289,
6998 3.2099992681699613513775259670214,
6999 3.24662674946606137764916854570219,
7000 3.28284687953866689817670991319787,
7001 3.31867291347259485044591136879087,
7002 3.35411740487202127264475726990106,
7003 3.38919225660177218727305224515862,
7004 3.42390876691942143189170489271753,
7005 3.45827767149820230182596660024454,
7006 3.49230918177808483937957161007792,
7007 3.5260130200285724149540352829756,
7008 3.55939845146044235497103883695448,
7009 3.59247431368364585025958062194665,
7010 3.62524904377393592090180712976368,
7011 3.65773070318071087226169680450936,
7012 3.68992700068237648299565823810245,
7013 3.72184531357268220291630708234186 };
7014 static doublereal xmax4[55] = { 1.1092649593311780079813740546678,
7015 1.05299572648705464724876659688996,
7016 1.0949715351434178709281698645813,
7017 1.15078388379719068145021100764647,
7018 1.2094863084718701596278219811869,
7019 1.26806623151369531323304177532868,
7020 1.32549784426476978866302826176202,
7021 1.38142537365039019558329304432581,
7022 1.43575531950773585146867625840552,
7023 1.48850442653629641402403231015299,
7024 1.53973611681876234549146350844736,
7025 1.58953193485272191557448229046492,
7026 1.63797820416306624705258190017418,
7027 1.68515974143594899185621942934906,
7028 1.73115699602477936547107755854868,
7029 1.77604489805513552087086912113251,
7030 1.81989256661534438347398400420601,
7031 1.86276344480103110090865609776681,
7032 1.90471563564740808542244678597105,
7033 1.94580231994751044968731427898046,
7034 1.98607219357764450634552790950067,
7035 2.02556989246317857340333585562678,
7036 2.06433638992049685189059517340452,
7037 2.10240936014742726236706004607473,
7038 2.13982350649113222745523925190532,
7039 2.17661085564771614285379929798896,
7040 2.21280102016879766322589373557048,
7041 2.2484214321456956597803794333791,
7042 2.28349755104077956674135810027654,
7043 2.31805304852593774867640120860446,
7044 2.35210997297725685169643559615022,
7045 2.38568889602346315560143377261814,
7046 2.41880904328694215730192284109322,
7047 2.45148841120796359750021227795539,
7048 2.48374387161372199992570528025315,
7049 2.5155912654873773953959098501893,
7050 2.54704548720896557684101746505398,
7051 2.57812056037881628390134077704127,
7052 2.60882970619319538196517982945269,
7053 2.63918540521920497868347679257107,
7054 2.66919945330942891495458446613851,
7055 2.69888301230439621709803756505788,
7056 2.72824665609081486737132853370048,
7057 2.75730041251405791603760003778285,
7058 2.78605380158311346185098508516203,
7059 2.81451587035387403267676338931454,
7060 2.84269522483114290814009184272637,
7061 2.87060005919012917988363332454033,
7062 2.89823818258367657739520912946934,
7063 2.92561704377132528239806135133273,
7064 2.95274375377994262301217318010209,
7065 2.97962510678256471794289060402033,
7066 3.00626759936182712291041810228171,
7067 3.03267744830655121818899164295959,
7068 3.05886060707437081434964933864149 };
7069 static doublereal xmax6[53] = { 1.21091229812484768570102219548814,
7070 1.11626917091567929907256116528817,
7071 1.1327140810290884106278510474203,
7072 1.1679452722668028753522098022171,
7073 1.20910611986279066645602153641334,
7074 1.25228283758701572089625983127043,
7075 1.29591971597287895911380446311508,
7076 1.3393138157481884258308028584917,
7077 1.3821288728999671920677617491385,
7078 1.42420414683357356104823573391816,
7079 1.46546895108549501306970087318319,
7080 1.50590085198398789708599726315869,
7081 1.54550385142820987194251585145013,
7082 1.58429644271680300005206185490937,
7083 1.62230484071440103826322971668038,
7084 1.65955905239130512405565733793667,
7085 1.69609056468292429853775667485212,
7086 1.73193098017228915881592458573809,
7087 1.7671112206990325429863426635397,
7088 1.80166107681586964987277458875667,
7089 1.83560897003644959204940535551721,
7090 1.86898184653271388435058371983316,
7091 1.90180515174518670797686768515502,
7092 1.93410285411785808749237200054739,
7093 1.96589749778987993293150856865539,
7094 1.99721027139062501070081653790635,
7095 2.02806108474738744005306947877164,
7096 2.05846864831762572089033752595401,
7097 2.08845055210580131460156962214748,
7098 2.11802334209486194329576724042253,
7099 2.14720259305166593214642386780469,
7100 2.17600297710595096918495785742803,
7101 2.20443832785205516555772788192013,
7102 2.2325216999457379530416998244706,
7103 2.2602654243075083168599953074345,
7104 2.28768115912702794202525264301585,
7105 2.3147799369092684021274946755348,
7106 2.34157220782483457076721300512406,
7107 2.36806787963276257263034969490066,
7108 2.39427635443992520016789041085844,
7109 2.42020656255081863955040620243062,
7110 2.44586699364757383088888037359254,
7111 2.47126572552427660024678584642791,
7112 2.49641045058324178349347438430311,
7113 2.52130850028451113942299097584818,
7114 2.54596686772399937214920135190177,
7115 2.5703922285006754089328998222275,
7116 2.59459096001908861492582631591134,
7117 2.61856915936049852435394597597773,
7118 2.64233265984385295286445444361827,
7119 2.66588704638685848486056711408168,
7120 2.68923766976735295746679957665724,
7121 2.71238965987606292679677228666411 };
7123 /* System generated locals */
7126 /* Local variables */
7127 static logical ldbg;
7128 static integer numax, ii;
7129 static doublereal bid;
7132 /* **********************************************************************
7137 /* Calculate the max of Jacobo polynoms multiplied by the weight on */
7138 /* (-1,1) for order 0,4,6 or Legendre. */
7142 /* LEGENDRE,APPROXIMATION,ERREUR. */
7144 /* INPUT ARGUMENTS : */
7145 /* ------------------ */
7146 /* NDGJAC: Nb of Jacobi coeff. of approximation. */
7147 /* IORDRE: Order of continuity (from -1 to 2) */
7149 /* OUTPUT ARGUMENTS : */
7150 /* ------------------- */
7151 /* XJACMX: Table of maximums of Jacobi polynoms. */
7153 /* COMMONS USED : */
7154 /* ---------------- */
7156 /* REFERENCES CALLED : */
7157 /* --------------------- */
7159 /* DESCRIPTION/NOTES/LIMITATIONS : */
7160 /* ----------------------------------- */
7163 /* ***********************************************************************
7165 /* Name of the routine */
7166 /* ----------------------------- Initialisations ------------------------
7169 ldbg = AdvApp2Var_SysBase::mnfndeb_() >= 3;
7171 AdvApp2Var_SysBase::mgenmsg_("MMA2JMX", 7L);
7174 numax = *ndgjac - ((*iordre + 1) << 1);
7175 if (*iordre == -1) {
7177 for (ii = 0; ii <= i__1; ++ii) {
7178 bid = (ii * 2. + 1.) / 2.;
7179 xjacmx[ii] = sqrt(bid);
7182 } else if (*iordre == 0) {
7184 for (ii = 0; ii <= i__1; ++ii) {
7185 xjacmx[ii] = xmax2[ii];
7188 } else if (*iordre == 1) {
7190 for (ii = 0; ii <= i__1; ++ii) {
7191 xjacmx[ii] = xmax4[ii];
7194 } else if (*iordre == 2) {
7196 for (ii = 0; ii <= i__1; ++ii) {
7197 xjacmx[ii] = xmax6[ii];
7202 /* ------------------------- The end ------------------------------------
7206 AdvApp2Var_SysBase::mgsomsg_("MMA2JMX", 7L);
7211 //=======================================================================
7212 //function : mma2moy_
7214 //=======================================================================
7215 int mma2moy_(integer *ndgumx,
7227 /* System generated locals */
7228 integer patjac_dim1, patjac_dim2, patjac_offset, i__1, i__2, i__3;
7230 /* Local variables */
7231 static logical ldbg;
7232 static integer minu, minv, idebu, idebv, ii, nd, jj;
7233 static doublereal bid0, bid1;
7236 /* **********************************************************************
7241 /* Calculate the average approximation error made when only */
7242 /* the coefficients of PATJAC of degree between */
7243 /* 2*(IORDRU+1) and MINDGU by U and 2*(IORDRV+1) and MINDGV by V are preserved. */
7247 /* LEGENDRE,APPROXIMATION, AVERAGE ERROR */
7249 /* INPUT ARGUMENTS : */
7250 /* ------------------ */
7251 /* NDGUMX: Dimension by U of table PATJAC. */
7252 /* NDGVMX: Dimension by V of table PATJAC. */
7253 /* NDIMEN: Dimension of the space. */
7254 /* MINDGU: Lower limit of the index by U of PATJAC coeff to be taken into account. */
7255 /* MAXDGU: Upper limit of the index by U of PATJAC coeff to be taken into account. */
7256 /* MINDGV: Lower limit of the index by V of PATJAC coeff to be taken into account. */
7257 /* MAXDGV: Upper limit of the index by V of PATJAC coeff to be taken into account. */
7258 /* IORDRU: Order of continuity by U provided by square PATJAC (from -1 to 2) */
7259 /* IORDRV: Order of continuity by V provided by square PATJAC (from -1 to 2) */
7260 /* PATJAC: Table of coeff. of the approximation square with */
7261 /* constraints of order IORDRU by U and IORDRV by V. */
7263 /* OUTPUT ARGUMENTS : */
7264 /* ------------------- */
7265 /* ERRMOY: Average error commited by preserving only the coeff of */
7266 /* PATJAC 2*(IORDRU+1) in MINDGU by U and 2*(IORDRV+1) in MINDGV by V. */
7268 /* COMMONS USED : */
7269 /* ---------------- */
7271 /* REFERENCES CALLED : */
7272 /* --------------------- */
7274 /* DESCRIPTION/NOTES/LIMITATIONS : */
7275 /* ----------------------------------- */
7276 /* Table PATJAC stores the coeff. Cij of */
7277 /* approximation square F(U,V). Indexes i and j show the degree by */
7278 /* U and by V of the base polynoms. These base polynoms are in the form: */
7280 /* ((1 - U*U)**(IORDRU+1)).J(i-2*(IORDRU+1)(U), where */
7282 /* polynom J(i-2*(IORDU+1)(U) is the Jacobi polynom of order */
7283 /* IORDRU+1 (the same by V by replacing U by V in the above expression). */
7285 /* The contribution to the average error of term Cij when */
7286 /* it is removed from PATJAC is Cij*Cij. */
7289 /* ***********************************************************************
7291 /* Name of the routine */
7294 /* ----------------------------- Initialisations ------------------------
7297 /* Parameter adjustments */
7298 patjac_dim1 = *ndgumx + 1;
7299 patjac_dim2 = *ndgvmx + 1;
7300 patjac_offset = patjac_dim1 * patjac_dim2;
7301 patjac -= patjac_offset;
7304 ldbg = AdvApp2Var_SysBase::mnfndeb_() >= 3;
7306 AdvApp2Var_SysBase::mgenmsg_("MMA2MOY", 7L);
7309 idebu = (*iordru + 1) << 1;
7310 idebv = (*iordrv + 1) << 1;
7311 minu = advapp_max(idebu,*mindgu);
7312 minv = advapp_max(idebv,*mindgv);
7316 /* ------------------ Calculation of the upper bound of the average error ------------ */
7317 /* -------------------- when the coeff. of indexes from MINDGU to MAXDGU ------ */
7318 /* ---------------- by U and of indexes from MINDGV to MAXDGV by V are removed -------------- */
7321 for (nd = 1; nd <= i__1; ++nd) {
7323 for (jj = minv; jj <= i__2; ++jj) {
7325 for (ii = idebu; ii <= i__3; ++ii) {
7326 bid1 = patjac[ii + (jj + nd * patjac_dim2) * patjac_dim1];
7327 bid0 += bid1 * bid1;
7336 for (nd = 1; nd <= i__1; ++nd) {
7338 for (jj = idebv; jj <= i__2; ++jj) {
7340 for (ii = minu; ii <= i__3; ++ii) {
7341 bid1 = patjac[ii + (jj + nd * patjac_dim2) * patjac_dim1];
7342 bid0 += bid1 * bid1;
7350 /* ----------------------- Calculation of the average error -------------
7354 *errmoy = sqrt(bid0);
7356 /* ------------------------- The end ------------------------------------
7360 AdvApp2Var_SysBase::mgsomsg_("MMA2MOY", 7L);
7365 //=======================================================================
7366 //function : mma2roo_
7368 //=======================================================================
7369 int AdvApp2Var_ApproxF2var::mma2roo_(integer *nbpntu,
7374 /* System generated locals */
7377 /* Local variables */
7378 static integer ii, ibb;
7380 /* **********************************************************************
7385 /* Return roots of Legendre for discretisations. */
7389 /* TOUS, AB_SPECIFI::CONTRAINTE&, DISCRETISATION, &POINT */
7391 /* INPUT ARGUMENTS : */
7392 /* ------------------ */
7393 /* NBPNTU: Nb of INTERNAL parameters of discretization BY U. */
7394 /* This is also the nb of root of the Legendre polynom where the discretization is done. */
7395 /* NBPNTV: Nb of INTERNAL parameters of discretization BY V. */
7396 /* This is also the nb of root of the Legendre polynom where the discretization is done. */
7398 /* OUTPUT ARGUMENTS : */
7399 /* ------------------- */
7400 /* UROOTL: Table of parameters of discretisation ON (-1,1) BY U.
7402 /* VROOTL: Table of parameters of discretisation ON (-1,1) BY V.
7405 /* COMMONS USED : */
7406 /* ---------------- */
7408 /* REFERENCES CALLED : */
7409 /* --------------------- */
7411 /* DESCRIPTION/NOTES/LIMITATIONS : */
7412 /* ----------------------------------- */
7415 /* **********************************************************************
7418 /* Name of the routine */
7421 /* Parameter adjustments */
7426 ibb = AdvApp2Var_SysBase::mnfndeb_();
7428 AdvApp2Var_SysBase::mgenmsg_("MMA2ROO", 7L);
7431 /* ---------------- Return the POSITIVE roots on U ------------------
7434 AdvApp2Var_MathBase::mmrtptt_(nbpntu, &urootl[(*nbpntu + 1) / 2 + 1]);
7436 for (ii = 1; ii <= i__1; ++ii) {
7437 urootl[ii] = -urootl[*nbpntu - ii + 1];
7440 if (*nbpntu % 2 == 1) {
7441 urootl[*nbpntu / 2 + 1] = 0.;
7444 /* ---------------- Return the POSITIVE roots on V ------------------
7447 AdvApp2Var_MathBase::mmrtptt_(nbpntv, &vrootl[(*nbpntv + 1) / 2 + 1]);
7449 for (ii = 1; ii <= i__1; ++ii) {
7450 vrootl[ii] = -vrootl[*nbpntv - ii + 1];
7453 if (*nbpntv % 2 == 1) {
7454 vrootl[*nbpntv / 2 + 1] = 0.;
7457 /* ------------------------------ The End -------------------------------
7461 AdvApp2Var_SysBase::mgsomsg_("MMA2ROO", 7L);
7465 //=======================================================================
7466 //function : mmmapcoe_
7468 //=======================================================================
7469 int mmmapcoe_(integer *ndim,
7479 /* System generated locals */
7480 integer somtab_dim1, somtab_offset, diftab_dim1, diftab_offset,
7481 crvjac_dim1, crvjac_offset, gsstab_dim1, i__1, i__2, i__3;
7483 /* Local variables */
7484 static integer igss, ikdeb;
7485 static doublereal bidon;
7486 static integer nd, ik, ir, nbroot, ibb;
7488 /* **********************************************************************
7493 /* Calculate the coefficients of polinomial approximation curve */
7494 /* of degree NDGJAC by the method of smallest squares starting from */
7495 /* the discretization of function on the roots of Legendre polynom */
7496 /* of degree NBPNTS. */
7500 /* FONCTION,APPROXIMATION,COEFFICIENT,POLYNOME */
7502 /* INPUT ARGUMENTS : */
7503 /* ------------------ */
7504 /* NDIM : Dimension of the space. */
7505 /* NDGJAC : Max Degree of the polynom of approximation. */
7506 /* The representation in the orthogonal base starts from degree */
7507 /* 0 to degree NDGJAC-2*(JORDRE+1). The polynomial base */
7508 /* is the base of Jacobi of order -1 (Legendre), 0, 1 and 2 */
7509 /* IORDRE : Order of the base of Jacobi (-1,0,1 or 2). Corresponds */
7510 /* to step of constraints, C0,C1 or C2. */
7511 /* NBPNTS : Degree of the polynom of Legendre on the roots which of */
7512 /* are calculated the coefficients of integration by */
7513 /* Gauss method. It is required to set NBPNTS=30,40,50 or 61 */
7514 /* and NDGJAC < NBPNTS. */
7515 /* SOMTAB : Table of F(ti)+F(-ti) with ti in ROOTAB. */
7516 /* DIFTAB : Table of F(ti)-F(-ti) with ti in ROOTAB. */
7517 /* GSSTAB(i,k) : Table of coefficients of integration by the Gauss method : */
7518 /* i varies from 0 to NBPNTS and */
7519 /* k varies from 0 to NDGJAC-2*(JORDRE+1). */
7521 /* OUTPUT ARGUMENTSE : */
7522 /* ------------------- */
7523 /* CRVJAC : Curve of approximation of FONCNP with eventually */
7524 /* taking into account of constraints at the extremities. */
7525 /* This curve is of degree NDGJAC. */
7527 /* COMMONS USED : */
7528 /* ---------------- */
7530 /* REFERENCES CALLED : */
7531 /* --------------------- */
7533 /* DESCRIPTION/NOTES/LIMITATIONS : */
7534 /* ------------------------------- */
7536 /* **********************************************************************
7539 /* Name of the routine */
7541 /* Parameter adjustments */
7542 crvjac_dim1 = *ndgjac + 1;
7543 crvjac_offset = crvjac_dim1;
7544 crvjac -= crvjac_offset;
7545 gsstab_dim1 = *nbpnts / 2 + 1;
7546 diftab_dim1 = *nbpnts / 2 + 1;
7547 diftab_offset = diftab_dim1;
7548 diftab -= diftab_offset;
7549 somtab_dim1 = *nbpnts / 2 + 1;
7550 somtab_offset = somtab_dim1;
7551 somtab -= somtab_offset;
7554 ibb = AdvApp2Var_SysBase::mnfndeb_();
7556 AdvApp2Var_SysBase::mgenmsg_("MMMAPCO", 7L);
7558 ikdeb = (*iordre + 1) << 1;
7559 nbroot = *nbpnts / 2;
7562 for (nd = 1; nd <= i__1; ++nd) {
7564 /* ----------------- Calculate the coefficients of even degree ----------
7568 for (ik = ikdeb; ik <= i__2; ik += 2) {
7572 for (ir = 1; ir <= i__3; ++ir) {
7573 bidon += somtab[ir + nd * somtab_dim1] * gsstab[ir + igss *
7577 crvjac[ik + nd * crvjac_dim1] = bidon;
7581 /* --------------- Calculate the coefficients of uneven degree ----------
7585 for (ik = ikdeb + 1; ik <= i__2; ik += 2) {
7589 for (ir = 1; ir <= i__3; ++ir) {
7590 bidon += diftab[ir + nd * diftab_dim1] * gsstab[ir + igss *
7594 crvjac[ik + nd * crvjac_dim1] = bidon;
7601 /* ------- Add terms connected to the supplementary root (0.D0) ------ */
7602 /* ----------- of Legendre polynom of uneven degree NBPNTS -----------
7605 if (*nbpnts % 2 == 0) {
7609 for (nd = 1; nd <= i__1; ++nd) {
7611 for (ik = ikdeb; ik <= i__2; ik += 2) {
7613 crvjac[ik + nd * crvjac_dim1] += somtab[nd * somtab_dim1] *
7614 gsstab[igss * gsstab_dim1];
7620 /* ------------------------------ The end -------------------------------
7625 AdvApp2Var_SysBase::mgsomsg_("MMMAPCO", 7L);
7629 //=======================================================================
7630 //function : mmaperm_
7632 //=======================================================================
7633 int mmaperm_(integer *ncofmx,
7641 /* System generated locals */
7642 integer crvjac_dim1, crvjac_offset, i__1, i__2;
7644 /* Local variables */
7645 static doublereal bidj;
7646 static integer i__, ia, nd, ncfcut, ibb;
7647 static doublereal bid;
7649 /* **********************************************************************
7654 /* Calculate the square root of the average quadratic error */
7655 /* of approximation done when only the */
7656 /* first NCFNEW coefficients of a curve of degree NCOEFF-1 */
7657 /* written in NORMALIZED Jacobi base of order 2*(IORDRE+1) are preserved. */
7661 /* LEGENDRE,POLYGONE,APPROXIMATION,ERREUR. */
7663 /* INPUT ARGUMENTS : */
7664 /* ------------------ */
7665 /* NCOFMX : Maximum degree of the curve. */
7666 /* NDIM : Dimension of the space. */
7667 /* NCOEFF : Degree +1 of the curve. */
7668 /* IORDRE : Order of constraint of continuity at the extremities. */
7669 /* CRVJAC : The curve the degree which of will be lowered. */
7670 /* NCFNEW : Degree +1 of the resulting polynom. */
7672 /* OUTPUT ARGUMENTS : */
7673 /* ------------------- */
7674 /* ERRMOY : Average precision of approximation. */
7676 /* COMMONS USED : */
7677 /* ---------------- */
7679 /* REFERENCES CALLED : */
7680 /* ----------------------- */
7682 /* DESCRIPTION/NOTES/LIMITATIONS : */
7683 /* ----------------------------------- */
7685 /* ***********************************************************************
7688 /* Name of the routine */
7690 /* Parameter adjustments */
7691 crvjac_dim1 = *ncofmx;
7692 crvjac_offset = crvjac_dim1 + 1;
7693 crvjac -= crvjac_offset;
7696 ibb = AdvApp2Var_SysBase::mnfndeb_();
7698 AdvApp2Var_SysBase::mgenmsg_("MMAPERM", 7L);
7701 /* --------- Minimum degree that can be reached : Stop at 1 or IA -------
7704 ia = (*iordre + 1) << 1;
7706 if (*ncfnew + 1 > ncfcut) {
7707 ncfcut = *ncfnew + 1;
7710 /* -------------- Elimination of coefficients of high degree ------------ */
7711 /* ----------- Loop on the series of Jacobi :NCFCUT --> NCOEFF --------- */
7716 for (nd = 1; nd <= i__1; ++nd) {
7718 for (i__ = ncfcut; i__ <= i__2; ++i__) {
7719 bidj = crvjac[i__ + nd * crvjac_dim1];
7726 /* ----------- Square Root of average quadratic error e -----------
7730 *errmoy = sqrt(bid);
7732 /* ------------------------------- The end ------------------------------
7736 AdvApp2Var_SysBase::mgsomsg_("MMAPERM", 7L);
7740 //=======================================================================
7741 //function : mmapptt_
7743 //=======================================================================
7744 int AdvApp2Var_ApproxF2var::mmapptt_(const integer *ndgjac,
7745 const integer *nbpnts,
7746 const integer *jordre,
7750 /* System generated locals */
7751 integer cgauss_dim1, i__1;
7753 /* Local variables */
7754 static integer kjac, iptt, ipdb0, infdg, iptdb, mxjac, ilong, ibb;
7756 /* **********************************************************************
7761 /* Load the elements required for integration by */
7762 /* Gauss method to obtain the coefficients in the base of
7763 /* Legendre of the approximation by the least squares of a */
7764 /* function. The elements are stored in commons MMAPGSS */
7765 /* (case without constraint), MMAPGS0 (constraints C0), MMAPGS1 */
7766 /* (constraints C1) and MMAPGS2 (constraints C2). */
7770 /* INTEGRATION,GAUSS,JACOBI */
7772 /* INPUT ARGUMENTS : */
7773 /* ------------------ */
7774 /* NDGJAC : Max degree of the polynom of approximation. */
7775 /* The representation in orthogonal base goes from degree
7776 /* 0 to degree NDGJAC-2*(JORDRE+1). The polynomial base */
7777 /* is the base of Jacobi of order -1 (Legendre), 0, 1 and 2 */
7778 /* NBPNTS : Degree of the polynom of Legendre on the roots which of */
7779 /* are calculated the coefficients of integration by the */
7780 /* method of Gauss. It is required that NBPNTS=8,10,15,20,25, */
7781 /* 30,40,50 or 61 and NDGJAC < NBPNTS. */
7782 /* JORDRE : Order of the base of Jacobi (-1,0,1 or 2). Corresponds */
7783 /* to step of constraints C0,C1 or C2. */
7785 /* OUTPUT ARGUMENTS : */
7786 /* ------------------- */
7787 /* CGAUSS(i,k) : Table of coefficients of integration by */
7788 /* Gauss method : i varies from 0 to the integer part */
7789 /* of NBPNTS/2 and k varies from 0 to NDGJAC-2*(JORDRE+1). */
7790 /* These are the coeff. of integration associated to */
7791 /* positive roots of the polynom of Legendre of degree */
7792 /* NBPNTS. CGAUSS(0,k) contains coeff. */
7793 /* of integration associated to root t = 0 when */
7794 /* NBPNTS is uneven. */
7795 /* IERCOD : Error code. */
7797 /* = 11 NBPNTS is not 8,10,15,20,25,30,40,50 or 61. */
7798 /* = 21 JORDRE is not -1,0,1 or 2. */
7799 /* = 31 NDGJAC is too great or too small. */
7801 /* COMMONS USED : */
7802 /* ---------------- */
7803 /* MMAPGSS,MMAPGS0,MMAPGS1,MMAPGS2. */
7804 /* ***********************************************************************
7806 /* Parameter adjustments */
7807 cgauss_dim1 = *nbpnts / 2 + 1;
7810 ibb = AdvApp2Var_SysBase::mnfndeb_();
7812 AdvApp2Var_SysBase::mgenmsg_("MMAPPTT", 7L);
7816 /* ------------------- Tests on the validity of inputs ----------------
7819 infdg = (*jordre + 1) << 1;
7820 if (*nbpnts != 8 && *nbpnts != 10 && *nbpnts != 15 && *nbpnts != 20 && *
7821 nbpnts != 25 && *nbpnts != 30 && *nbpnts != 40 && *nbpnts != 50 &&
7826 if (*jordre < -1 || *jordre > 2) {
7830 if (*ndgjac >= *nbpnts || *ndgjac < infdg) {
7834 /* --------------- Calculation of the start pointer following NBPNTS -----------
7839 iptdb += (8 - infdg) << 2;
7842 iptdb += (10 - infdg) * 5;
7845 iptdb += (15 - infdg) * 7;
7848 iptdb += (20 - infdg) * 10;
7851 iptdb += (25 - infdg) * 12;
7854 iptdb += (30 - infdg) * 15;
7857 iptdb += (40 - infdg) * 20;
7860 iptdb += (50 - infdg) * 25;
7865 ipdb0 = ipdb0 + (14 - infdg) / 2 + 1;
7868 ipdb0 = ipdb0 + (24 - infdg) / 2 + 1;
7871 /* ------------------ Choice of the common depending on JORDRE -------------
7874 if (*jordre == -1) {
7887 /* ---------------- Common MMAPGSS (case without constraints) ----------------
7891 ilong = *nbpnts / 2 << 3;
7893 for (kjac = 0; kjac <= i__1; ++kjac) {
7894 iptt = iptdb + kjac * (*nbpnts / 2) + 1;
7895 AdvApp2Var_SysBase::mcrfill_(&ilong,
7896 (char *)&mmapgss_.gslxjs[iptt - 1],
7897 (char *)&cgauss[kjac * cgauss_dim1 + 1]);
7900 /* --> Case when the number of points is uneven. */
7901 if (*nbpnts % 2 == 1) {
7904 for (kjac = 0; kjac <= i__1; kjac += 2) {
7905 cgauss[kjac * cgauss_dim1] = mmapgss_.gsl0js[iptt - 1];
7910 for (kjac = 1; kjac <= i__1; kjac += 2) {
7911 cgauss[kjac * cgauss_dim1] = 0.;
7917 /* ---------------- Common MMAPGS0 (case with constraints C0) -------------
7921 mxjac = *ndgjac - infdg;
7922 ilong = *nbpnts / 2 << 3;
7924 for (kjac = 0; kjac <= i__1; ++kjac) {
7925 iptt = iptdb + kjac * (*nbpnts / 2) + 1;
7926 AdvApp2Var_SysBase::mcrfill_(&ilong,
7927 (char *)&mmapgs0_.gslxj0[iptt - 1],
7928 (char *)&cgauss[kjac * cgauss_dim1 + 1]);
7931 /* --> Case when the number of points is uneven. */
7932 if (*nbpnts % 2 == 1) {
7935 for (kjac = 0; kjac <= i__1; kjac += 2) {
7936 cgauss[kjac * cgauss_dim1] = mmapgs0_.gsl0j0[iptt - 1];
7941 for (kjac = 1; kjac <= i__1; kjac += 2) {
7942 cgauss[kjac * cgauss_dim1] = 0.;
7948 /* ---------------- Common MMAPGS1 (case with constraints C1) -------------
7952 mxjac = *ndgjac - infdg;
7953 ilong = *nbpnts / 2 << 3;
7955 for (kjac = 0; kjac <= i__1; ++kjac) {
7956 iptt = iptdb + kjac * (*nbpnts / 2) + 1;
7957 AdvApp2Var_SysBase::mcrfill_(&ilong,
7958 (char *)&mmapgs1_.gslxj1[iptt - 1],
7959 (char *)&cgauss[kjac * cgauss_dim1 + 1]);
7962 /* --> Case when the number of points is uneven. */
7963 if (*nbpnts % 2 == 1) {
7966 for (kjac = 0; kjac <= i__1; kjac += 2) {
7967 cgauss[kjac * cgauss_dim1] = mmapgs1_.gsl0j1[iptt - 1];
7972 for (kjac = 1; kjac <= i__1; kjac += 2) {
7973 cgauss[kjac * cgauss_dim1] = 0.;
7979 /* ---------------- Common MMAPGS2 (case with constraints C2) -------------
7983 mxjac = *ndgjac - infdg;
7984 ilong = *nbpnts / 2 << 3;
7986 for (kjac = 0; kjac <= i__1; ++kjac) {
7987 iptt = iptdb + kjac * (*nbpnts / 2) + 1;
7988 AdvApp2Var_SysBase::mcrfill_(&ilong,
7989 (char *)&mmapgs2_.gslxj2[iptt - 1],
7990 (char *)&cgauss[kjac * cgauss_dim1 + 1]);
7993 /* --> Cas of uneven number of points. */
7994 if (*nbpnts % 2 == 1) {
7997 for (kjac = 0; kjac <= i__1; kjac += 2) {
7998 cgauss[kjac * cgauss_dim1] = mmapgs2_.gsl0j2[iptt - 1];
8003 for (kjac = 1; kjac <= i__1; kjac += 2) {
8004 cgauss[kjac * cgauss_dim1] = 0.;
8010 /* ------------------------- Return the error code --------------
8012 /* --> NBPNTS is not OK */
8016 /* --> JORDRE is not OK */
8020 /* --> NDGJAC is not OK */
8025 /* -------------------------------- The end -----------------------------
8030 AdvApp2Var_SysBase::maermsg_("MMAPPTT", iercod, 7L);
8033 AdvApp2Var_SysBase::mgsomsg_("MMAPPTT", 7L);
8039 //=======================================================================
8040 //function : mmjacpt_
8042 //=======================================================================
8043 int mmjacpt_(const integer *ndimen,
8044 const integer *ncoefu,
8045 const integer *ncoefv,
8046 const integer *iordru,
8047 const integer *iordrv,
8048 const doublereal *ptclgd,
8052 /* System generated locals */
8053 integer ptccan_dim1, ptccan_dim2, ptccan_offset, ptclgd_dim1, ptclgd_dim2,
8054 ptclgd_offset, ptcaux_dim1, ptcaux_dim2, ptcaux_dim3,
8055 ptcaux_offset, i__1, i__2, i__3;
8057 /* Local variables */
8058 static integer kdim, nd, ii, jj, ibb;
8060 /* ***********************************************************************
8065 /* Passage from canonical to Jacobi base for a */
8066 /* "square" in a space of arbitrary dimension. */
8070 /* SMOOTHING,BASE,LEGENDRE */
8073 /* INPUT ARGUMENTS : */
8074 /* ------------------ */
8075 /* NDIMEN : Dimension of the space. */
8076 /* NCOEFU : Degree+1 by U. */
8077 /* NCOEFV : Degree+1 by V. */
8078 /* IORDRU : Order of Jacobi polynoms by U. */
8079 /* IORDRV : Order of Jacobi polynoms by V. */
8080 /* PTCLGD : The square in the Jacobi base. */
8082 /* OUTPUT ARGUMENTS : */
8083 /* ------------------- */
8084 /* PTCAUX : Auxilliary space. */
8085 /* PTCCAN : The square in the canonic base (-1,1) */
8087 /* COMMONS USED : */
8088 /* ---------------- */
8090 /* APPLIED REFERENCES : */
8091 /* ----------------------- */
8093 /* DESCRIPTION/NOTES/LIMITATIONS : */
8094 /* ----------------------------------- */
8095 /* Cancels and replaces MJACPC */
8097 /* *********************************************************************
8099 /* Name of the routine */
8102 /* Parameter adjustments */
8103 ptccan_dim1 = *ncoefu;
8104 ptccan_dim2 = *ncoefv;
8105 ptccan_offset = ptccan_dim1 * (ptccan_dim2 + 1) + 1;
8106 ptccan -= ptccan_offset;
8107 ptcaux_dim1 = *ncoefv;
8108 ptcaux_dim2 = *ncoefu;
8109 ptcaux_dim3 = *ndimen;
8110 ptcaux_offset = ptcaux_dim1 * (ptcaux_dim2 * (ptcaux_dim3 + 1) + 1) + 1;
8111 ptcaux -= ptcaux_offset;
8112 ptclgd_dim1 = *ncoefu;
8113 ptclgd_dim2 = *ncoefv;
8114 ptclgd_offset = ptclgd_dim1 * (ptclgd_dim2 + 1) + 1;
8115 ptclgd -= ptclgd_offset;
8118 ibb = AdvApp2Var_SysBase::mnfndeb_();
8120 AdvApp2Var_SysBase::mgenmsg_("MMJACPT", 7L);
8123 /* Passage into canonical by u. */
8125 kdim = *ndimen * *ncoefv;
8126 AdvApp2Var_MathBase::mmjaccv_((integer *)ncoefu,
8129 (doublereal *)&ptclgd[ptclgd_offset],
8130 (doublereal *)&ptcaux[ptcaux_offset],
8131 (doublereal *)&ptccan[ptccan_offset]);
8133 /* Swapping of u and v. */
8136 for (nd = 1; nd <= i__1; ++nd) {
8138 for (jj = 1; jj <= i__2; ++jj) {
8140 for (ii = 1; ii <= i__3; ++ii) {
8141 ptcaux[jj + (ii + (nd + ptcaux_dim3) * ptcaux_dim2) *
8142 ptcaux_dim1] = ptccan[ii + (jj + nd * ptccan_dim2) *
8151 /* Passage into canonical by v. */
8153 kdim = *ndimen * *ncoefu;
8154 AdvApp2Var_MathBase::mmjaccv_((integer *)ncoefv,
8157 (doublereal *)&ptcaux[((ptcaux_dim3 + 1) * ptcaux_dim2 + 1) * ptcaux_dim1 + 1],
8158 (doublereal *)&ptccan[ptccan_offset],
8159 (doublereal *)&ptcaux[(((ptcaux_dim3 << 1) + 1) * ptcaux_dim2 + 1) * ptcaux_dim1 + 1]);
8161 /* Swapping of u and v. */
8164 for (nd = 1; nd <= i__1; ++nd) {
8166 for (jj = 1; jj <= i__2; ++jj) {
8168 for (ii = 1; ii <= i__3; ++ii) {
8169 ptccan[ii + (jj + nd * ptccan_dim2) * ptccan_dim1] = ptcaux[
8170 jj + (ii + (nd + (ptcaux_dim3 << 1)) * ptcaux_dim2) *
8179 /* ---------------------------- THAT'S ALL FOLKS ------------------------
8183 AdvApp2Var_SysBase::mgsomsg_("MMJACPT", 7L);