1 // Copyright (c) 1999-2012 OPEN CASCADE SAS
3 // The content of this file is subject to the Open CASCADE Technology Public
4 // License Version 6.5 (the "License"). You may not use the content of this file
5 // except in compliance with the License. Please obtain a copy of the License
6 // at http://www.opencascade.org and read it completely before using this file.
8 // The Initial Developer of the Original Code is Open CASCADE S.A.S., having its
9 // main offices at: 1, place des Freres Montgolfier, 78280 Guyancourt, France.
11 // The Original Code and all software distributed under the License is
12 // distributed on an "AS IS" basis, without warranty of any kind, and the
13 // Initial Developer hereby disclaims all such warranties, including without
14 // limitation, any warranties of merchantability, fitness for a particular
15 // purpose or non-infringement. Please see the License for the specific terms
16 // and conditions governing the rights and limitations under the License.
18 // AdvApp2Var_ApproxF2var.cxx
20 #include <AdvApp2Var_SysBase.hxx>
21 #include <AdvApp2Var_MathBase.hxx>
22 #include <AdvApp2Var_Data_f2c.hxx>
23 #include <AdvApp2Var_Data.hxx>
24 #include <AdvApp2Var_ApproxF2var.hxx>
28 int mmjacpt_(const integer *ndimen,
29 const integer *ncoefu,
30 const integer *ncoefv,
31 const integer *iordru,
32 const integer *iordrv,
33 const doublereal *ptclgd,
40 int mma2ce2_(integer *numdec,
75 int mma2cfu_(integer *ndujac,
87 int mma2cfv_(integer *ndvjac,
97 int mma2er1_(integer *ndjacu,
113 int mma2er2_(integer *ndjacu,
132 int mma2moy_(integer *ndgumx,
145 int mma2ds2_(integer *ndimen,
148 const AdvApp2Var_EvaluatorFunc2Var& foncnp,
166 int mma1fdi_(integer *ndimen,
168 const AdvApp2Var_EvaluatorFunc2Var& foncnp,
183 int mma1cdi_(integer *ndimen,
195 int mma1jak_(integer *ndimen,
205 int mma1cnt_(integer *ndimen,
214 int mma1fer_(integer *ndimen,
229 int mma1noc_(doublereal *dfuvin,
240 int mmmapcoe_(integer *ndim,
250 int mmaperm_(integer *ncofmx,
259 #define mmapgss_1 mmapgss_
260 #define mmapgs0_1 mmapgs0_
261 #define mmapgs1_1 mmapgs1_
262 #define mmapgs2_1 mmapgs2_
264 //=======================================================================
265 //function : mma1cdi_
267 //=======================================================================
268 int mma1cdi_(integer *ndimen,
282 /* System generated locals */
283 integer contr1_dim1, contr1_offset, contr2_dim1, contr2_offset,
284 somtab_dim1, somtab_offset, diftab_dim1, diftab_offset,
285 fpntab_dim1, fpntab_offset, hermit_dim1, hermit_offset, i__1,
288 /* Local variables */
289 integer nroo2, ncfhe, nd, ii, kk;
290 integer ibb, kkm, kkp;
291 doublereal bid1, bid2, bid3 = 0.;
293 /* **********************************************************************
297 /* Discretisation on the parameters of interpolation polynomes */
298 /* constraints of order IORDRE. */
302 /* ALL, AB_SPECIFI::CONTRAINTE&, DISCRETISATION, &POINT */
304 /* INPUT ARGUMENTS : */
305 /* ------------------ */
306 /* NDIMEN: Space dimension. */
307 /* NBROOT: Number of INTERNAL discretisation parameters. */
308 /* It is also the root number Legendre polynome where */
309 /* the discretization is performed. */
310 /* ROOTLG: Table of discretization parameters ON (-1,1). */
311 /* IORDRE: Order of constraint imposed to the extremities of the iso. */
312 /* = 0, the extremities of the iso are calculated */
313 /* = 1, additionally, the 1st derivative in the direction */
314 /* of the iso is calculated. */
315 /* = 2, additionally, the 2nd derivative in the direction */
316 /* of the iso is calculated. */
317 /* CONTR1: Contains, if IORDRE>=0, values IORDRE+1 in TTABLE(0)
319 /* (1st extremity) of derivatives of F(Uc,Ve) or F(Ue,Vc), */
321 /* CONTR2: Contains, if IORDRE>=0, values IORDRE+1 in */
322 /* TTABLE(NBROOT+1) (2nd extremity) of: */
323 /* If ISOFAV=1, derived of order IDERIV by U, derived */
324 /* ordre 0 to IORDRE by V of F(Uc,Ve) or Uc=TCONST */
325 /* (fixed iso value) and Ve is the fixed extremity. */
326 /* If ISOFAV=2, derivative of order IDERIV by V, derivative */
327 /* of order 0 to IORDRE by U of F(Ue,Vc) or Vc=TCONST */
328 /* (fixed iso value) and Ue is the fixed extremity. */
330 /* SOMTAB: Table of NBROOT/2 sums of 2 index points */
331 /* NBROOT-II+1 and II, for II = 1, NBROOT/2. */
332 /* DIFTAB: Table of NBROOT/2 differences of 2 index points */
333 /* NBROOT-II+1 and II, for II = 1, NBROOT/2. */
335 /* OUTPUT ARGUMENTS : */
336 /* ------------------- */
337 /* SOMTAB: Table of NBROOT/2 sums of 2 index points */
338 /* NBROOT-II+1 and II, for II = 1, NBROOT/2 */
339 /* DIFTAB: Table of NBROOT/2 differences of 2 index points */
340 /* NBROOT-II+1 and II, for II = 1, NBROOT/2 */
341 /* FPNTAB: Auxiliary table. */
342 /* HERMIT: Table of coeff. 2*(IORDRE+1) Hermite polynoms */
343 /* of degree 2*IORDRE+1. */
344 /* IERCOD: Error code, */
345 /* = 0, Everythig is OK */
346 /* = 1, The value of IORDRE is out of (0,2) */
348 /* ---------------- */
350 /* REFERENCES CALLED : */
351 /* ----------------------- */
353 /* DESCRIPTION/NOTES/LIMITATIONS : */
354 /* ----------------------------------- */
355 /* The results of discretization are arranged in 2 tables */
356 /* SOMTAB and DIFTAB to earn time during the */
357 /* calculation of coefficients of the approximation curve. */
359 /* If NBROOT is uneven in SOMTAB(0,*) and DIFTAB(0,*) one stores */
360 /* the values of the median root of Legendre (0.D0 in (-1,1)). */
362 /* **********************************************************************
365 /* Name of the routine */
368 /* Parameter adjustments */
369 diftab_dim1 = *nbroot / 2 + 1;
370 diftab_offset = diftab_dim1;
371 diftab -= diftab_offset;
372 somtab_dim1 = *nbroot / 2 + 1;
373 somtab_offset = somtab_dim1;
374 somtab -= somtab_offset;
376 hermit_dim1 = (*iordre << 1) + 2;
377 hermit_offset = hermit_dim1;
378 hermit -= hermit_offset;
379 fpntab_dim1 = *nbroot;
380 fpntab_offset = fpntab_dim1 + 1;
381 fpntab -= fpntab_offset;
382 contr2_dim1 = *ndimen;
383 contr2_offset = contr2_dim1 + 1;
384 contr2 -= contr2_offset;
385 contr1_dim1 = *ndimen;
386 contr1_offset = contr1_dim1 + 1;
387 contr1 -= contr1_offset;
390 ibb = AdvApp2Var_SysBase::mnfndeb_();
392 AdvApp2Var_SysBase::mgenmsg_("MMA1CDI", 7L);
396 /* --- Recuperate 2*(IORDRE+1) coeff of 2*(IORDRE+1) of Hermite polynom ---
399 AdvApp2Var_ApproxF2var::mma1her_(iordre, &hermit[hermit_offset], iercod);
404 /* ------------------- Discretization of Hermite polynoms -----------
407 ncfhe = (*iordre + 1) << 1;
409 for (ii = 1; ii <= i__1; ++ii) {
411 for (kk = 1; kk <= i__2; ++kk) {
412 AdvApp2Var_MathBase::mmmpocur_(&ncfhe, &c__1, &ncfhe, &hermit[ii * hermit_dim1], &
413 rootlg[kk], &fpntab[kk + ii * fpntab_dim1]);
419 /* ---- Discretizations of boundary polynoms are taken ----
424 for (nd = 1; nd <= i__1; ++nd) {
426 for (ii = 1; ii <= i__2; ++ii) {
427 bid1 = contr1[nd + ii * contr1_dim1];
428 bid2 = contr2[nd + ii * contr2_dim1];
430 for (kk = 1; kk <= i__3; ++kk) {
431 kkm = nroo2 - kk + 1;
432 bid3 = bid1 * fpntab[kkm + ((ii << 1) - 1) * fpntab_dim1] +
433 bid2 * fpntab[kkm + (ii << 1) * fpntab_dim1];
434 somtab[kk + nd * somtab_dim1] -= bid3;
435 diftab[kk + nd * diftab_dim1] += bid3;
439 for (kk = 1; kk <= i__3; ++kk) {
440 kkp = (*nbroot + 1) / 2 + kk;
441 bid3 = bid1 * fpntab[kkp + ((ii << 1) - 1) * fpntab_dim1] +
442 bid2 * fpntab[kkp + (ii << 1) * fpntab_dim1];
443 somtab[kk + nd * somtab_dim1] -= bid3;
444 diftab[kk + nd * diftab_dim1] -= bid3;
452 /* ------------ Cas when discretization is done on the roots of a -----------
454 /* ---------- Legendre polynom of uneven degree, 0 is root --------
457 if (*nbroot % 2 == 1) {
459 for (nd = 1; nd <= i__1; ++nd) {
461 for (ii = 1; ii <= i__2; ++ii) {
462 bid3 = fpntab[nroo2 + 1 + ((ii << 1) - 1) * fpntab_dim1] *
463 contr1[nd + ii * contr1_dim1] + fpntab[nroo2 + 1 + (
464 ii << 1) * fpntab_dim1] * contr2[nd + ii *
468 somtab[nd * somtab_dim1] -= bid3;
469 diftab[nd * diftab_dim1] -= bid3;
476 /* ------------------------------ The End -------------------------------
478 /* --> IORDRE is not in the authorized zone. */
485 AdvApp2Var_SysBase::mgsomsg_("MMA1CDI", 7L);
490 //=======================================================================
491 //function : mma1cnt_
493 //=======================================================================
494 int mma1cnt_(integer *ndimen,
502 /* System generated locals */
503 integer contr1_dim1, contr1_offset, contr2_dim1, contr2_offset,
504 hermit_dim1, hermit_offset, crvjac_dim1, crvjac_offset, i__1,
507 /* Local variables */
508 integer nd, ii, jj, ibb;
512 /* ***********************************************************************
517 /* Add constraint to polynom. */
521 /* ALL,AB_SPECIFI::COURE&,APPROXIMATION,ADDITION,&CONSTRAINT */
523 /* INPUT ARGUMENTS : */
524 /* -------------------- */
525 /* NDIMEN: Dimension of the space */
526 /* IORDRE: Order of constraint. */
527 /* CONTR1: pt of constraint in -1, from order 0 to IORDRE. */
528 /* CONTR2: Pt of constraint in +1, from order 0 to IORDRE. */
529 /* HERMIT: Table of Hermit polynoms of order IORDRE. */
530 /* CRVJAV: Curve of approximation in Jacobi base. */
532 /* OUTPUT ARGUMENTS : */
533 /* --------------------- */
534 /* CRVJAV: Curve of approximation in Jacobi base */
535 /* to which the polynom of interpolation of constraints is added. */
538 /* ------------------ */
541 /* REFERENCES CALLED : */
542 /* --------------------- */
545 /* DESCRIPTION/NOTES/LIMITATIONS : */
546 /* ----------------------------------- */
549 /* ***********************************************************************
552 /* ***********************************************************************
554 /* Name of the routine */
556 /* ***********************************************************************
558 /* INITIALISATIONS */
559 /* ***********************************************************************
562 /* Parameter adjustments */
563 hermit_dim1 = (*iordre << 1) + 2;
564 hermit_offset = hermit_dim1;
565 hermit -= hermit_offset;
566 contr2_dim1 = *ndimen;
567 contr2_offset = contr2_dim1 + 1;
568 contr2 -= contr2_offset;
569 contr1_dim1 = *ndimen;
570 contr1_offset = contr1_dim1 + 1;
571 contr1 -= contr1_offset;
572 crvjac_dim1 = *ndgjac + 1;
573 crvjac_offset = crvjac_dim1;
574 crvjac -= crvjac_offset;
577 ibb = AdvApp2Var_SysBase::mnfndeb_();
579 AdvApp2Var_SysBase::mgenmsg_("MMA1CNT", 7L);
582 /* ***********************************************************************
585 /* ***********************************************************************
589 for (nd = 1; nd <= i__1; ++nd) {
590 i__2 = (*iordre << 1) + 1;
591 for (ii = 0; ii <= i__2; ++ii) {
594 for (jj = 1; jj <= i__3; ++jj) {
595 bid = bid + contr1[nd + jj * contr1_dim1] *
596 hermit[ii + ((jj << 1) - 1) * hermit_dim1] +
597 contr2[nd + jj * contr2_dim1] * hermit[ii + (jj << 1) * hermit_dim1];
600 crvjac[ii + nd * crvjac_dim1] = bid;
606 /* ***********************************************************************
608 /* RETURN CALLING PROGRAM */
609 /* ***********************************************************************
613 AdvApp2Var_SysBase::mgsomsg_("MMA1CNT", 7L);
619 //=======================================================================
620 //function : mma1fdi_
622 //=======================================================================
623 int mma1fdi_(integer *ndimen,
625 const AdvApp2Var_EvaluatorFunc2Var& foncnp,
639 /* System generated locals */
640 integer fpntab_dim1, somtab_dim1, somtab_offset, diftab_dim1,
641 diftab_offset, contr1_dim1, contr1_offset, contr2_dim1,
642 contr2_offset, i__1, i__2;
645 /* Local variables */
646 integer ideb, ifin, nroo2, ideru, iderv;
648 integer ii, nd, ibb, iim, nbp, iip;
649 doublereal bid1, bid2;
651 /* **********************************************************************
656 /* DiscretiZation of a non-polynomial function F(U,V) or of */
657 /* its derivative with fixed isoparameter. */
661 /* ALL, AB_SPECIFI::FONCTION&, DISCRETISATION, &POINT */
663 /* INPUT ARGUMENTS : */
664 /* ------------------ */
665 /* NDIMEN: Space dimension. */
666 /* UVFONC: Limits of the path of definition by U and by V of the approximated function */
667 /* FONCNP: The NAME of the non-polynomial function to be approximated */
668 /* (external program). */
669 /* ISOFAV: Fixed isoparameter for the discretization; */
670 /* = 1, discretization with fixed U and variable V. */
671 /* = 2, discretization with fixed V and variable U. */
672 /* TCONST: Iso value is also fixed. */
673 /* NBROOT: Number of INTERNAL discretization parameters. */
674 /* (if there are constraints, 2 extremities should be added).
676 /* This is also the root number of the Legendre polynom where */
677 /* the discretization is done. */
678 /* TTABLE: Table of discretization parameters and of 2 extremities */
679 /* (Respectively (-1, NBROOT Legendre roots,1) */
680 /* reframed within the adequate interval. */
681 /* IORDRE: Order of constraint imposed on the extremities of the iso. */
682 /* (If Iso-U, it is necessary to calculate the derivatives by V and vice */
684 /* = 0, the extremities of the iso are calculated. */
685 /* = 1, additionally the 1st derivative in the direction of the iso is calculated */
686 /* = 2, additionally the 2nd derivative in the direction of the iso is calculated */
687 /* IDERIV: Order of derivative transversal to fixed iso (If Iso-U=Uc */
688 /* is fixed, the derivative of order IDERIV is discretized by U of */
689 /* F(Uc,v). Same if iso-V is fixed). */
690 /* Varies from 0 (positioning) to 2 (2nd derivative). */
692 /* OUTPUT ARGUMENTS : */
693 /* ------------------- */
694 /* FPNTAB: Auxiliary table.
695 SOMTAB: Table of NBROOT/2 sums of 2 index points */
696 /* NBROOT-II+1 and II, for II = 1, NBROOT/2 */
697 /* DIFTAB: Table of NBROOT/2 differences of 2 index points */
698 /* NBROOT-II+1 and II, for II = 1, NBROOT/2 */
699 /* CONTR1: Contains, if IORDRE>=0, values IORDRE+1 in TTABLE(0)
701 /* (1st extremity) of derivatives of F(Uc,Ve) or F(Ue,Vc), */
703 /* CONTR2: Contains, if IORDRE>=0, values IORDRE+1 in */
704 /* TTABLE(NBROOT+1) (2nd extremity) of: */
705 /* If ISOFAV=1, derived of order IDERIV by U, derived */
706 /* ordre 0 to IORDRE by V of F(Uc,Ve) or Uc=TCONST */
707 /* (fixed iso value) and Ve is the fixed extremity. */
708 /* If ISOFAV=2, derivative of order IDERIV by V, derivative */
709 /* of order 0 to IORDRE by U of F(Ue,Vc) or Vc=TCONST */
710 /* (fixed iso value) and Ue is the fixed extremity. */
711 /* IERCOD: Error code > 100; Pb in evaluation of FONCNP, */
712 /* the returned error code is equal to error code of FONCNP + 100. */
715 /* ---------------- */
717 /* REFERENCES CALLED : */
718 /* ----------------------- */
720 /* DESCRIPTION/NOTES/LIMITATIONS : */
721 /* ----------------------------------- */
722 /* The results of discretization are arranged in 2 tables */
723 /* SOMTAB and DIFTAB to earn time during the */
724 /* calculation of coefficients of the approximation curve. */
726 /* If NBROOT is uneven in SOMTAB(0,*) and DIFTAB(0,*) one stores */
727 /* the values of the median root of Legendre (0.D0 in (-1,1)). */
729 /* Function F(u,v) defined in UVFONC is reparameterized in */
730 /* (-1,1)x(-1,1). Then 1st and 2nd derivatives are renormalized. */
733 /* **********************************************************************
736 /* Name of the routine */
739 /* Parameter adjustments */
741 diftab_dim1 = *nbroot / 2 + 1;
742 diftab_offset = diftab_dim1;
743 diftab -= diftab_offset;
744 somtab_dim1 = *nbroot / 2 + 1;
745 somtab_offset = somtab_dim1;
746 somtab -= somtab_offset;
747 fpntab_dim1 = *ndimen;
749 contr2_dim1 = *ndimen;
750 contr2_offset = contr2_dim1 + 1;
751 contr2 -= contr2_offset;
752 contr1_dim1 = *ndimen;
753 contr1_offset = contr1_dim1 + 1;
754 contr1 -= contr1_offset;
757 ibb = AdvApp2Var_SysBase::mnfndeb_();
759 AdvApp2Var_SysBase::mgenmsg_("MMA1FDI", 7L);
763 /* --------------- Definition of the nb of points to calculate --------------
765 /* --> If constraints, the limits are also taken */
769 /* --> Otherwise, only Legendre roots (reframed) are used
775 /* --> Nb of point to calculate. */
776 nbp = ifin - ideb + 1;
779 /* --------------- Determination of the order of global derivation --------
781 /* --> ISOFAV takes only values 1 or 2. */
782 /* if Iso-U, derive by U of order IDERIV */
786 d__1 = (uvfonc[4] - uvfonc[3]) / 2.;
787 renor = AdvApp2Var_MathBase::pow__di(&d__1, ideriv);
788 /* if Iso-V, derive by V of order IDERIV */
792 d__1 = (uvfonc[6] - uvfonc[5]) / 2.;
793 renor = AdvApp2Var_MathBase::pow__di(&d__1, ideriv);
796 /* ----------- Discretization on roots of the ---------------
798 /* ---------------------- Legendre polynom of degree NBROOT -------------------
801 (*const_cast <AdvApp2Var_EvaluatorFunc2Var*> (&foncnp)).Evaluate (ndimen,
810 &fpntab[ideb * fpntab_dim1 + 1],
816 for (nd = 1; nd <= i__1; ++nd) {
818 for (ii = 1; ii <= i__2; ++ii) {
819 iip = (*nbroot + 1) / 2 + ii;
820 iim = nroo2 - ii + 1;
821 bid1 = fpntab[nd + iim * fpntab_dim1];
822 bid2 = fpntab[nd + iip * fpntab_dim1];
823 somtab[ii + nd * somtab_dim1] = renor * (bid2 + bid1);
824 diftab[ii + nd * diftab_dim1] = renor * (bid2 - bid1);
830 /* ------------ Case when discretisation is done on roots of a ----
832 /* ---------- Legendre polynom of uneven degree, 0 is root --------
835 if (*nbroot % 2 == 1) {
837 for (nd = 1; nd <= i__1; ++nd) {
838 somtab[nd * somtab_dim1] = renor * fpntab[nd + (nroo2 + 1) *
840 diftab[nd * diftab_dim1] = renor * fpntab[nd + (nroo2 + 1) *
846 for (nd = 1; nd <= i__1; ++nd) {
847 somtab[nd * somtab_dim1] = 0.;
848 diftab[nd * diftab_dim1] = 0.;
853 /* --------------------- Take into account constraints ----------------
857 /* --> Recover already calculated extremities. */
859 for (nd = 1; nd <= i__1; ++nd) {
860 contr1[nd + contr1_dim1] = renor * fpntab[nd];
861 contr2[nd + contr2_dim1] = renor * fpntab[nd + (*nbroot + 1) *
865 /* --> Nb of points to calculate/call to FONCNP */
867 /* If Iso-U, derive by V till order IORDRE */
869 /* --> Factor of normalisation 1st derivative. */
870 bid1 = (uvfonc[6] - uvfonc[5]) / 2.;
872 for (iderv = 1; iderv <= i__1; ++iderv) {
873 (*const_cast <AdvApp2Var_EvaluatorFunc2Var*> (&foncnp)).Evaluate (
874 ndimen, &uvfonc[3], &uvfonc[5], isofav, tconst, &
875 nbp, ttable, &ideru, &iderv, &contr1[(iderv + 1) *
876 contr1_dim1 + 1], iercod);
883 for (iderv = 1; iderv <= i__1; ++iderv) {
884 (*const_cast <AdvApp2Var_EvaluatorFunc2Var*> (&foncnp)).Evaluate (
885 ndimen, &uvfonc[3], &uvfonc[5], isofav, tconst, &
886 nbp, &ttable[*nbroot + 1], &ideru, &iderv, &contr2[(
887 iderv + 1) * contr2_dim1 + 1], iercod);
893 /* If Iso-V, derive by U till order IORDRE */
895 /* --> Factor of normalization 1st derivative. */
896 bid1 = (uvfonc[4] - uvfonc[3]) / 2.;
898 for (ideru = 1; ideru <= i__1; ++ideru) {
899 (*const_cast <AdvApp2Var_EvaluatorFunc2Var*> (&foncnp)).Evaluate (
900 ndimen, &uvfonc[3], &uvfonc[5], isofav, tconst, &
901 nbp, ttable, &ideru, &iderv, &contr1[(ideru + 1) *
902 contr1_dim1 + 1], iercod);
909 for (ideru = 1; ideru <= i__1; ++ideru) {
910 (*const_cast <AdvApp2Var_EvaluatorFunc2Var*> (&foncnp)).Evaluate (
911 ndimen, &uvfonc[3], &uvfonc[5], isofav, tconst, &
912 nbp, &ttable[*nbroot + 1], &ideru, &iderv, &contr2[(
913 ideru + 1) * contr2_dim1 + 1], iercod);
921 /* ------------------------- Normalization of derivatives -------------
923 /* (The function is redefined on (-1,1)*(-1,1)) */
926 for (ii = 1; ii <= i__1; ++ii) {
929 for (nd = 1; nd <= i__2; ++nd) {
930 contr1[nd + (ii + 1) * contr1_dim1] *= bid2;
931 contr2[nd + (ii + 1) * contr2_dim1] *= bid2;
938 /* ------------------------------ The end -------------------------------
944 AdvApp2Var_SysBase::maermsg_("MMA1FDI", iercod, 7L);
947 AdvApp2Var_SysBase::mgsomsg_("MMA1FDI", 7L);
952 //=======================================================================
953 //function : mma1fer_
955 //=======================================================================
956 int mma1fer_(integer *,//ndimen,
970 /* System generated locals */
971 integer crvjac_dim1, crvjac_offset, i__1, i__2;
973 /* Local variables */
974 integer idim, ncfja, ncfnw, ndses, ii, kk, ibb, ier;
978 /* ***********************************************************************
983 /* Calculate the degree and the errors of approximation of a border. */
987 /* TOUS,AB_SPECIFI :: COURBE&,TRONCATURE, &PRECISION */
989 /* INPUT ARGUMENTS : */
990 /* -------------------- */
992 /* NDIMEN: Total Dimension of the space (sum of dimensions of sub-spaces) */
993 /* NBSESP: Number of "independent" sub-spaces. */
994 /* NDIMSE: Table of dimensions of sub-spaces. */
995 /* IORDRE: Order of constraint at the extremities of the border */
996 /* -1 = no constraints, */
997 /* 0 = constraints of passage to limits (i.e. C0), */
998 /* 1 = C0 + constraintes of 1st derivatives (i.e. C1), */
999 /* 2 = C1 + constraintes of 2nd derivatives (i.e. C2). */
1000 /* NDGJAC: Degree of development in series to use for the calculation
1001 /* in the base of Jacobi. */
1002 /* CRVJAC: Table of coeff. of the curve of approximation in the */
1003 /* base of Jacobi. */
1004 /* NCFLIM: Max number of coeff of the polynomial curve */
1005 /* of approximation (should be above or equal to */
1006 /* 2*IORDRE+2 and below or equal to 50). */
1007 /* EPSAPR: Table of errors of approximations that cannot be passed, */
1008 /* sub-space by sub-space. */
1010 /* OUTPUT ARGUMENTS : */
1011 /* --------------------- */
1012 /* YCVMAX: Auxiliary Table. */
1013 /* ERRMAX: Table of errors (sub-space by sub-space) */
1014 /* MAXIMUM made in the approximation of FONCNP by */
1016 /* ERRMOY: Table of errors (sub-space by sub-space) */
1017 /* AVERAGE made in the approximation of FONCNP by */
1019 /* NCOEFF: Number of significative coeffs. of the calculated "curve". */
1020 /* IERCOD: Error code */
1022 /* =-1, warning, required tolerance can't be */
1023 /* met with coefficients NFCLIM. */
1024 /* = 1, order of constraints (IORDRE) is not within authorised values */
1027 /* COMMONS USED : */
1028 /* ------------------ */
1030 /* REFERENCES CALLED : */
1031 /* --------------------- */
1033 /* DESCRIPTION/NOTES/LIMITATIONS : */
1034 /* ----------------------------------- */
1036 /* **********************************************************************
1039 /* Name of the routine */
1042 /* Parameter adjustments */
1048 crvjac_dim1 = *ndgjac + 1;
1049 crvjac_offset = crvjac_dim1;
1050 crvjac -= crvjac_offset;
1053 ibb = AdvApp2Var_SysBase::mnfndeb_();
1055 AdvApp2Var_SysBase::mgenmsg_("MMA1FER", 7L);
1060 ncfja = *ndgjac + 1;
1062 /* ------------ Calculate the degree of the curve and of the Max error --------
1064 /* -------------- of approximation for all sub-spaces --------
1068 for (ii = 1; ii <= i__1; ++ii) {
1071 /* ------------ cutting of coeff. and calculation of Max error -------
1074 AdvApp2Var_MathBase::mmtrpjj_(&ncfja, &ndses, &ncfja, &epsapr[ii], iordre, &crvjac[idim *
1075 crvjac_dim1], &ycvmax[1], &errmax[ii], &ncfnw);
1077 /* ******************************************************************
1079 /* ------------- If precision OK, calculate the average error -------
1081 /* ******************************************************************
1084 if (ncfnw <= *ncflim) {
1085 mmaperm_(&ncfja, &ndses, &ncfja, iordre, &crvjac[idim *
1086 crvjac_dim1], &ncfnw, &errmoy[ii]);
1087 *ncoeff = advapp_max(ncfnw,*ncoeff);
1089 /* ------------- Set the declined coefficients to 0.D0 -----------
1092 nbr0 = *ncflim - ncfnw;
1095 for (kk = 1; kk <= i__2; ++kk) {
1096 AdvApp2Var_SysBase::mvriraz_(&nbr0,
1097 &crvjac[ncfnw + (idim + kk - 1) * crvjac_dim1]);
1103 /* **************************************************************
1105 /* ------------------- If required precision can't be reached----
1107 /* **************************************************************
1112 /* ------------------------- calculate the Max error ------------
1115 AdvApp2Var_MathBase::mmaperx_(&ncfja, &ndses, &ncfja, iordre, &crvjac[idim *
1116 crvjac_dim1], ncflim, &ycvmax[1], &errmax[ii], &ier);
1121 /* -------------------- nb of coeff to be returned -------------
1126 /* ------------------- and calculation of the average error ----
1129 mmaperm_(&ncfja, &ndses, &ncfja, iordre, &crvjac[idim *
1130 crvjac_dim1], ncflim, &errmoy[ii]);
1138 /* ------------------------------ The end -------------------------------
1140 /* --> The order of constraints is not within autorized values. */
1147 AdvApp2Var_SysBase::maermsg_("MMA1FER", iercod, 7L);
1150 AdvApp2Var_SysBase::mgsomsg_("MMA1FER", 7L);
1156 //=======================================================================
1157 //function : mma1her_
1159 //=======================================================================
1160 int AdvApp2Var_ApproxF2var::mma1her_(const integer *iordre,
1164 /* System generated locals */
1165 integer hermit_dim1, hermit_offset;
1167 /* Local variables */
1172 /* **********************************************************************
1177 /* Calculate 2*(IORDRE+1) Hermit polynoms of degree 2*IORDRE+1 */
1182 /* ALL, AB_SPECIFI::CONTRAINTE&, INTERPOLATION, &POLYNOME */
1184 /* INPUT ARGUMENTS : */
1185 /* ------------------ */
1186 /* IORDRE: Order of constraint. */
1187 /* = 0, Polynom of interpolation of order C0 on (-1,1). */
1188 /* = 1, Polynom of interpolation of order C0 and C1 on (-1,1). */
1189 /* = 2, Polynom of interpolation of order C0, C1 and C2 on (-1,1).
1192 /* OUTPUT ARGUMENTS : */
1193 /* ------------------- */
1194 /* HERMIT: Table of 2*IORDRE+2 coeff. of each of 2*(IORDRE+1) */
1195 /* HERMIT polynom. */
1196 /* IERCOD: Error code, */
1198 /* = 1, required order of constraint is not managed here. */
1199 /* COMMONS USED : */
1200 /* ---------------- */
1202 /* REFERENCES CALLED : */
1203 /* ----------------------- */
1205 /* DESCRIPTION/NOTES/LIMITATIONS : */
1206 /* ----------------------------------- */
1207 /* The part of HERMIT(*,2*i+j) table where j=1 or 2 and i=0 to IORDRE,
1208 /* contains the coefficients of the polynom of degree 2*IORDRE+1 */
1209 /* such as ALL values in -1 and in +1 of this polynom and its */
1210 /* derivatives till order of derivation IORDRE are NULL, */
1211 /* EXCEPT for the derivative of order i: */
1212 /* - valued 1 in -1 if j=1 */
1213 /* - valued 1 in +1 if j=2. */
1215 /* **********************************************************************
1218 /* Name of the routine */
1221 /* Parameter adjustments */
1222 hermit_dim1 = (*iordre + 1) << 1;
1223 hermit_offset = hermit_dim1 + 1;
1224 hermit -= hermit_offset;
1227 ibb = AdvApp2Var_SysBase::mnfndeb_();
1229 AdvApp2Var_SysBase::mgenmsg_("MMA1HER", 7L);
1233 /* --- Recover (IORDRE+2) coeff of 2*(IORDRE+1) Hermit polynoms --
1237 hermit[hermit_dim1 + 1] = .5;
1238 hermit[hermit_dim1 + 2] = -.5;
1240 hermit[(hermit_dim1 << 1) + 1] = .5;
1241 hermit[(hermit_dim1 << 1) + 2] = .5;
1242 } else if (*iordre == 1) {
1243 hermit[hermit_dim1 + 1] = .5;
1244 hermit[hermit_dim1 + 2] = -.75;
1245 hermit[hermit_dim1 + 3] = 0.;
1246 hermit[hermit_dim1 + 4] = .25;
1248 hermit[(hermit_dim1 << 1) + 1] = .5;
1249 hermit[(hermit_dim1 << 1) + 2] = .75;
1250 hermit[(hermit_dim1 << 1) + 3] = 0.;
1251 hermit[(hermit_dim1 << 1) + 4] = -.25;
1253 hermit[hermit_dim1 * 3 + 1] = .25;
1254 hermit[hermit_dim1 * 3 + 2] = -.25;
1255 hermit[hermit_dim1 * 3 + 3] = -.25;
1256 hermit[hermit_dim1 * 3 + 4] = .25;
1258 hermit[(hermit_dim1 << 2) + 1] = -.25;
1259 hermit[(hermit_dim1 << 2) + 2] = -.25;
1260 hermit[(hermit_dim1 << 2) + 3] = .25;
1261 hermit[(hermit_dim1 << 2) + 4] = .25;
1262 } else if (*iordre == 2) {
1263 hermit[hermit_dim1 + 1] = .5;
1264 hermit[hermit_dim1 + 2] = -.9375;
1265 hermit[hermit_dim1 + 3] = 0.;
1266 hermit[hermit_dim1 + 4] = .625;
1267 hermit[hermit_dim1 + 5] = 0.;
1268 hermit[hermit_dim1 + 6] = -.1875;
1270 hermit[(hermit_dim1 << 1) + 1] = .5;
1271 hermit[(hermit_dim1 << 1) + 2] = .9375;
1272 hermit[(hermit_dim1 << 1) + 3] = 0.;
1273 hermit[(hermit_dim1 << 1) + 4] = -.625;
1274 hermit[(hermit_dim1 << 1) + 5] = 0.;
1275 hermit[(hermit_dim1 << 1) + 6] = .1875;
1277 hermit[hermit_dim1 * 3 + 1] = .3125;
1278 hermit[hermit_dim1 * 3 + 2] = -.4375;
1279 hermit[hermit_dim1 * 3 + 3] = -.375;
1280 hermit[hermit_dim1 * 3 + 4] = .625;
1281 hermit[hermit_dim1 * 3 + 5] = .0625;
1282 hermit[hermit_dim1 * 3 + 6] = -.1875;
1284 hermit[(hermit_dim1 << 2) + 1] = -.3125;
1285 hermit[(hermit_dim1 << 2) + 2] = -.4375;
1286 hermit[(hermit_dim1 << 2) + 3] = .375;
1287 hermit[(hermit_dim1 << 2) + 4] = .625;
1288 hermit[(hermit_dim1 << 2) + 5] = -.0625;
1289 hermit[(hermit_dim1 << 2) + 6] = -.1875;
1291 hermit[hermit_dim1 * 5 + 1] = .0625;
1292 hermit[hermit_dim1 * 5 + 2] = -.0625;
1293 hermit[hermit_dim1 * 5 + 3] = -.125;
1294 hermit[hermit_dim1 * 5 + 4] = .125;
1295 hermit[hermit_dim1 * 5 + 5] = .0625;
1296 hermit[hermit_dim1 * 5 + 6] = -.0625;
1298 hermit[hermit_dim1 * 6 + 1] = .0625;
1299 hermit[hermit_dim1 * 6 + 2] = .0625;
1300 hermit[hermit_dim1 * 6 + 3] = -.125;
1301 hermit[hermit_dim1 * 6 + 4] = -.125;
1302 hermit[hermit_dim1 * 6 + 5] = .0625;
1303 hermit[hermit_dim1 * 6 + 6] = .0625;
1308 /* ------------------------------ The End -------------------------------
1311 AdvApp2Var_SysBase::maermsg_("MMA1HER", iercod, 7L);
1313 AdvApp2Var_SysBase::mgsomsg_("MMA1HER", 7L);
1317 //=======================================================================
1318 //function : mma1jak_
1320 //=======================================================================
1321 int mma1jak_(integer *ndimen,
1331 /* System generated locals */
1332 integer somtab_dim1, somtab_offset, diftab_dim1, diftab_offset,
1333 crvjac_dim1, crvjac_offset, cgauss_dim1;
1335 /* Local variables */
1338 /* **********************************************************************
1343 /* Calculate the curve of approximation of a non-polynomial function */
1344 /* in the base of Jacobi. */
1348 /* FUNCTION,DISCRETISATION,APPROXIMATION,CONSTRAINT,CURVE,JACOBI */
1350 /* INPUT ARGUMENTS : */
1351 /* ------------------ */
1352 /* NDIMEN: Total dimension of the space (sum of dimensions */
1353 /* of sub-spaces) */
1354 /* NBROOT: Nb of points of discretization of the iso, extremities not
1356 /* IORDRE: Order of constraint at the extremities of the boundary */
1357 /* -1 = no constraints, */
1358 /* 0 = constraints of passage of limits (i.e. C0), */
1359 /* 1 = C0 + constraints of 1st derivatives (i.e. C1), */
1360 /* 2 = C1 + constraints of 2nd derivatives (i.e. C2). */
1361 /* NDGJAC: Degree of development in series to be used for calculation in the
1362 /* base of Jacobi. */
1364 /* OUTPUT ARGUMENTS : */
1365 /* ------------------- */
1366 /* CRVJAC : Curve of approximation of FONCNP with (eventually) */
1367 /* taking into account of constraints at the extremities. */
1368 /* This curve is of degree NDGJAC. */
1369 /* IERCOD : Error code : */
1370 /* 0 = All is ok. */
1371 /* 33 = Pb to return data of du block data */
1372 /* of coeff. of integration by GAUSS method. */
1373 /* by program MMAPPTT. */
1375 /* COMMONS USED : */
1376 /* ---------------- */
1378 /* REFERENCES CALLED : */
1379 /* ----------------------- */
1381 /* **********************************************************************
1384 /* Name of the routine */
1386 /* Parameter adjustments */
1387 diftab_dim1 = *nbroot / 2 + 1;
1388 diftab_offset = diftab_dim1;
1389 diftab -= diftab_offset;
1390 somtab_dim1 = *nbroot / 2 + 1;
1391 somtab_offset = somtab_dim1;
1392 somtab -= somtab_offset;
1393 crvjac_dim1 = *ndgjac + 1;
1394 crvjac_offset = crvjac_dim1;
1395 crvjac -= crvjac_offset;
1396 cgauss_dim1 = *nbroot / 2 + 1;
1399 ibb = AdvApp2Var_SysBase::mnfndeb_();
1401 AdvApp2Var_SysBase::mgenmsg_("MMA1JAK", 7L);
1405 /* ----------------- Recover coeffs of integration by Gauss -----------
1408 AdvApp2Var_ApproxF2var::mmapptt_(ndgjac, nbroot, iordre, cgauss, iercod);
1414 /* --------------- Calculate the curve in the base of Jacobi -----------
1417 mmmapcoe_(ndimen, ndgjac, iordre, nbroot, &somtab[somtab_offset], &diftab[
1418 diftab_offset], cgauss, &crvjac[crvjac_offset]);
1420 /* ------------------------------ The End -------------------------------
1425 AdvApp2Var_SysBase::maermsg_("MMA1JAK", iercod, 7L);
1428 AdvApp2Var_SysBase::mgsomsg_("MMA1JAK", 7L);
1433 //=======================================================================
1434 //function : mma1noc_
1436 //=======================================================================
1437 int mma1noc_(doublereal *dfuvin,
1446 /* System generated locals */
1450 /* Local variables */
1451 doublereal rider, riord;
1454 /* **********************************************************************
1459 /* Normalization of constraints of derivatives, defined on DFUVIN */
1460 /* on block DUVOUT. */
1464 /* ALL, AB_SPECIFI::VECTEUR&,DERIVEE&,NORMALISATION,&VECTEUR */
1466 /* INPUT ARGUMENTS : */
1467 /* ------------------ */
1468 /* DFUVIN: Limits of the block of definition by U and by V where
1470 /* constraints CNTRIN are defined. */
1471 /* NDIMEN: Dimension of the space. */
1472 /* IORDRE: Order of constraint imposed at the extremities of the iso. */
1473 /* (if Iso-U, it is necessary to calculate derivatives by V and vice */
1475 /* = 0, the extremities of the iso are calculated */
1476 /* = 1, additionally the 1st derivative in the direction */
1477 /* of the iso is calculated */
1478 /* = 2, additionally the 2nd derivative in the direction */
1479 /* of the iso is calculated */
1480 /* CNTRIN: Contains, if IORDRE>=0, IORDRE+1 derivatives */
1481 /* of order IORDRE of F(Uc,v) or of F(u,Vc), following the */
1482 /* value of ISOFAV, RENORMALIZED by u and v in (-1,1). */
1483 /* DUVOUT: Limits of the block of definition by U and by V where the */
1484 /* constraints CNTOUT will be defined. */
1485 /* ISOFAV: Isoparameter fixed for the discretization; */
1486 /* = 1, discretization with fixed U=Uc and variable V. */
1487 /* = 2, discretization with fixed V=Vc and variable U. */
1488 /* IDERIV: Ordre de derivee transverse a l'iso fixee (Si Iso-U=Uc */
1489 /* is fixed, the derivative of order IDERIV is discretized by U */
1490 /* of F(Uc,v). The same if iso-V is fixed). */
1491 /* Varies from (positioning) to 2 (2nd derivative). */
1493 /* OUTPUT ARGUMENTS : */
1494 /* ------------------- */
1495 /* CNTOUT: Contains, if IORDRE>=0, IORDRE+1 derivatives */
1496 /* of order IORDRE of F(Uc,v) or of F(u,Vc), depending on the */
1497 /* value of ISOFAV, RENORMALIZED for u and v in DUVOUT. */
1499 /* COMMONS USED : */
1500 /* ---------------- */
1502 /* REFERENCES CALLED : */
1503 /* --------------------- */
1505 /* DESCRIPTION/NOTES/LIMITATIONS : */
1506 /* ------------------------------- */
1507 /* CNTRIN can be an output/input argument, */
1510 /* CALL MMA1NOC(DFUVIN,NDIMEN,IORDRE,CNTRIN,DUVOUT */
1511 /* 1 ,ISOFAV,IDERIV,CNTRIN) */
1515 /* **********************************************************************
1518 /* Name of the routine */
1521 /* Parameter adjustments */
1528 ibb = AdvApp2Var_SysBase::mnfndeb_();
1530 AdvApp2Var_SysBase::mgenmsg_("MMA1NOC", 7L);
1533 /* --------------- Determination of coefficients of normalization -------
1537 d__1 = (dfuvin[4] - dfuvin[3]) / (duvout[4] - duvout[3]);
1538 rider = AdvApp2Var_MathBase::pow__di(&d__1, ideriv);
1539 d__1 = (dfuvin[6] - dfuvin[5]) / (duvout[6] - duvout[5]);
1540 riord = AdvApp2Var_MathBase::pow__di(&d__1, iordre);
1543 d__1 = (dfuvin[6] - dfuvin[5]) / (duvout[6] - duvout[5]);
1544 rider = AdvApp2Var_MathBase::pow__di(&d__1, ideriv);
1545 d__1 = (dfuvin[4] - dfuvin[3]) / (duvout[4] - duvout[3]);
1546 riord = AdvApp2Var_MathBase::pow__di(&d__1, iordre);
1549 /* ------------- Renormalization of the vector of constraint ---------------
1552 bid = rider * riord;
1554 for (nd = 1; nd <= i__1; ++nd) {
1555 cntout[nd] = bid * cntrin[nd];
1559 /* ------------------------------ The end -------------------------------
1563 AdvApp2Var_SysBase::mgsomsg_("MMA1NOC", 7L);
1568 //=======================================================================
1569 //function : mma1nop_
1571 //=======================================================================
1572 int mma1nop_(integer *nbroot,
1580 /* System generated locals */
1583 /* Local variables */
1584 doublereal alinu, blinu, alinv, blinv;
1587 /* ***********************************************************************
1592 /* Normalization of parameters of an iso, starting from */
1593 /* parametric block and parameters on (-1,1). */
1597 /* TOUS,AB_SPECIFI :: ISO&,POINT&,NORMALISATION,&POINT,&ISO */
1599 /* INPUT ARGUMENTS : */
1600 /* -------------------- */
1601 /* NBROOT: Nb of points of discretisation INSIDE the iso */
1602 /* defined on (-1,1). */
1603 /* ROOTLG: Table of discretization parameters on )-1,1( */
1605 /* UVFONC: Block of definition of the iso */
1606 /* ISOFAV: = 1, this is iso-u; =2, this is iso-v. */
1608 /* OUTPUT ARGUMENTS : */
1609 /* --------------------- */
1610 /* TTABLE: Table of parameters renormalized on UVFONC of the iso.
1612 /* IERCOD: = 0, OK */
1613 /* = 1, ISOFAV is out of allowed values. */
1616 /* **********************************************************************
1618 /* Name of the routine */
1621 /* Parameter adjustments */
1626 ibb = AdvApp2Var_SysBase::mnfndeb_();
1628 AdvApp2Var_SysBase::mgenmsg_("MMA1NOP", 7L);
1631 alinu = (uvfonc[4] - uvfonc[3]) / 2.;
1632 blinu = (uvfonc[4] + uvfonc[3]) / 2.;
1633 alinv = (uvfonc[6] - uvfonc[5]) / 2.;
1634 blinv = (uvfonc[6] + uvfonc[5]) / 2.;
1637 ttable[0] = uvfonc[5];
1639 for (ii = 1; ii <= i__1; ++ii) {
1640 ttable[ii] = alinv * rootlg[ii] + blinv;
1643 ttable[*nbroot + 1] = uvfonc[6];
1644 } else if (*isofav == 2) {
1645 ttable[0] = uvfonc[3];
1647 for (ii = 1; ii <= i__1; ++ii) {
1648 ttable[ii] = alinu * rootlg[ii] + blinu;
1651 ttable[*nbroot + 1] = uvfonc[4];
1658 /* ------------------------------ THE END -------------------------------
1667 AdvApp2Var_SysBase::maermsg_("MMA1NOP", iercod, 7L);
1670 AdvApp2Var_SysBase::mgsomsg_("MMA1NOP", 7L);
1677 //=======================================================================
1678 //function : mma2ac1_
1680 //=======================================================================
1681 int AdvApp2Var_ApproxF2var::mma2ac1_(integer const *ndimen,
1682 integer const *mxujac,
1683 integer const *mxvjac,
1684 integer const *iordru,
1685 integer const *iordrv,
1686 doublereal const *contr1,
1687 doublereal const * contr2,
1688 doublereal const *contr3,
1689 doublereal const *contr4,
1690 doublereal const *uhermt,
1691 doublereal const *vhermt,
1695 /* System generated locals */
1696 integer contr1_dim1, contr1_dim2, contr1_offset, contr2_dim1, contr2_dim2,
1697 contr2_offset, contr3_dim1, contr3_dim2, contr3_offset,
1698 contr4_dim1, contr4_dim2, contr4_offset, uhermt_dim1,
1699 uhermt_offset, vhermt_dim1, vhermt_offset, patjac_dim1,
1700 patjac_dim2, patjac_offset, i__1, i__2, i__3, i__4, i__5;
1702 /* Local variables */
1705 doublereal bidu1, bidu2, bidv1, bidv2;
1706 integer ioru1, iorv1, ii, nd, jj, ku, kv;
1707 doublereal cnt1, cnt2, cnt3, cnt4;
1709 /* **********************************************************************
1714 /* Add polynoms of edge constraints. */
1718 /* TOUS,AB_SPECIFI::POINT&,CONTRAINTE&,ADDITION,&POLYNOME */
1720 /* INPUT ARGUMENTS : */
1721 /* ------------------ */
1722 /* NDIMEN: Dimension of the space. */
1723 /* MXUJAC: Max degree of the polynom of approximation by U. The */
1724 /* representation in the orthogonal base starts from degree */
1725 /* 0 to degree MXUJAC-2*(IORDRU+1). The polynomial base is the */
1726 /* base of Jacobi of order -1 (Legendre), 0, 1 or 2 */
1727 /* MXVJAC: Max degree of the polynom of approximation by V. The */
1728 /* representation in the orthogonal base starts from degree */
1729 /* 0 to degree MXUJAC-2*(IORDRU+1). The polynomial base is the */
1730 /* base of Jacobi of order -1 (Legendre), 0, 1 or 2 */
1731 /* IORDRU: Order of the base of Jacobi (-1,0,1 or 2) by U. Corresponds */
1732 /* to the step of constraints: C0, C1 or C2. */
1733 /* IORDRV: Order of the base of Jacobi (-1,0,1 or 2) by V. Corresponds */
1734 /* to the step of constraints: C0, C1 or C2. */
1735 /* CONTR1: Contains, if IORDRU and IORDRV>=0, the values at the */
1736 /* extremities of F(U0,V0) and its derivatives. */
1737 /* CONTR2: Contains, if IORDRU and IORDRV>=0, the values at the */
1738 /* extremities of F(U1,V0) and its derivatives. */
1739 /* CONTR3: Contains, if IORDRU and IORDRV>=0, the values at the */
1740 /* extremities of F(U0,V1) and its derivatives. */
1741 /* CONTR4: Contains, if IORDRU and IORDRV>=0, the values at the */
1742 /* extremities of F(U1,V1) and its derivatives. */
1743 /* UHERMT: Coeff. of Hermit polynoms of order IORDRU. */
1744 /* VHERMT: Coeff. of Hermit polynoms of order IORDRV. */
1745 /* PATJAC: Table of coefficients of the polynom P(u,v) of approximation */
1746 /* of F(u,v) WITHOUT taking into account the constraints. */
1748 /* OUTPUT ARGUMENTS : */
1749 /* ------------------- */
1750 /* PATJAC: Table of coefficients of the polynom P(u,v) by approximation */
1751 /* of F(u,v) WITH taking into account of constraints. */
1753 /* **********************************************************************
1755 /* Name of the routine */
1757 /* --------------------------- Initialization --------------------------
1760 /* Parameter adjustments */
1761 patjac_dim1 = *mxujac + 1;
1762 patjac_dim2 = *mxvjac + 1;
1763 patjac_offset = patjac_dim1 * patjac_dim2;
1764 patjac -= patjac_offset;
1765 uhermt_dim1 = (*iordru << 1) + 2;
1766 uhermt_offset = uhermt_dim1;
1767 uhermt -= uhermt_offset;
1768 vhermt_dim1 = (*iordrv << 1) + 2;
1769 vhermt_offset = vhermt_dim1;
1770 vhermt -= vhermt_offset;
1771 contr4_dim1 = *ndimen;
1772 contr4_dim2 = *iordru + 2;
1773 contr4_offset = contr4_dim1 * (contr4_dim2 + 1) + 1;
1774 contr4 -= contr4_offset;
1775 contr3_dim1 = *ndimen;
1776 contr3_dim2 = *iordru + 2;
1777 contr3_offset = contr3_dim1 * (contr3_dim2 + 1) + 1;
1778 contr3 -= contr3_offset;
1779 contr2_dim1 = *ndimen;
1780 contr2_dim2 = *iordru + 2;
1781 contr2_offset = contr2_dim1 * (contr2_dim2 + 1) + 1;
1782 contr2 -= contr2_offset;
1783 contr1_dim1 = *ndimen;
1784 contr1_dim2 = *iordru + 2;
1785 contr1_offset = contr1_dim1 * (contr1_dim2 + 1) + 1;
1786 contr1 -= contr1_offset;
1789 ldbg = AdvApp2Var_SysBase::mnfndeb_() >= 3;
1791 AdvApp2Var_SysBase::mgenmsg_("MMA2AC1", 7L);
1794 /* ------------ SUBTRACTION OF ANGULAR CONSTRAINTS -------------------
1797 ioru1 = *iordru + 1;
1798 iorv1 = *iordrv + 1;
1799 ndgu = (*iordru << 1) + 1;
1800 ndgv = (*iordrv << 1) + 1;
1803 for (jj = 1; jj <= i__1; ++jj) {
1805 for (ii = 1; ii <= i__2; ++ii) {
1807 for (nd = 1; nd <= i__3; ++nd) {
1808 cnt1 = contr1[nd + (ii + jj * contr1_dim2) * contr1_dim1];
1809 cnt2 = contr2[nd + (ii + jj * contr2_dim2) * contr2_dim1];
1810 cnt3 = contr3[nd + (ii + jj * contr3_dim2) * contr3_dim1];
1811 cnt4 = contr4[nd + (ii + jj * contr4_dim2) * contr4_dim1];
1813 for (kv = 0; kv <= i__4; ++kv) {
1814 bidv1 = vhermt[kv + ((jj << 1) - 1) * vhermt_dim1];
1815 bidv2 = vhermt[kv + (jj << 1) * vhermt_dim1];
1817 for (ku = 0; ku <= i__5; ++ku) {
1818 bidu1 = uhermt[ku + ((ii << 1) - 1) * uhermt_dim1];
1819 bidu2 = uhermt[ku + (ii << 1) * uhermt_dim1];
1820 patjac[ku + (kv + nd * patjac_dim2) * patjac_dim1] =
1821 patjac[ku + (kv + nd * patjac_dim2) *
1822 patjac_dim1] - bidu1 * bidv1 * cnt1 - bidu2 *
1823 bidv1 * cnt2 - bidu1 * bidv2 * cnt3 - bidu2 *
1836 /* ------------------------------ The end -------------------------------
1840 AdvApp2Var_SysBase::mgsomsg_("MMA2AC1", 7L);
1845 //=======================================================================
1846 //function : mma2ac2_
1848 //=======================================================================
1849 int AdvApp2Var_ApproxF2var::mma2ac2_(const integer *ndimen,
1850 const integer *mxujac,
1851 const integer *mxvjac,
1852 const integer *iordrv,
1853 const integer *nclimu,
1854 const integer *ncfiv1,
1855 const doublereal *crbiv1,
1856 const integer *ncfiv2,
1857 const doublereal *crbiv2,
1858 const doublereal *vhermt,
1862 /* System generated locals */
1863 integer crbiv1_dim1, crbiv1_dim2, crbiv1_offset, crbiv2_dim1, crbiv2_dim2,
1864 crbiv2_offset, patjac_dim1, patjac_dim2, patjac_offset,
1865 vhermt_dim1, vhermt_offset, i__1, i__2, i__3, i__4;
1867 /* Local variables */
1869 integer ndgv1, ndgv2, ii, jj, nd, kk;
1870 doublereal bid1, bid2;
1872 /* **********************************************************************
1877 /* Add polynoms of constraints */
1881 /* FUNCTION,APPROXIMATION,COEFFICIENT,POLYNOM */
1883 /* INPUT ARGUMENTS : */
1884 /* ------------------ */
1885 /* NDIMEN: Dimension of the space. */
1886 /* MXUJAC: Max degree of the polynom of approximation by U. The */
1887 /* representation in the orthogonal base starts from degree */
1888 /* 0 to degree MXUJAC-2*(IORDRU+1). The polynomial base is the */
1889 /* base of Jacobi of order -1 (Legendre), 0, 1 or 2 */
1890 /* MXVJAC: Max degree of the polynom of approximation by V. The */
1891 /* representation in the orthogonal base starts from degree */
1892 /* 0 to degree MXUJAC-2*(IORDRU+1). The polynomial base is the */
1893 /* base of Jacobi of order -1 (Legendre), 0, 1 or 2 */
1894 /* IORDRV: Order of the base of Jacobi (-1,0,1 or 2) by V. Corresponds */
1895 /* to the step of constraints: C0, C1 or C2. */
1896 /* NCLIMU LIMIT nb of coeff by u of the solution P(u,v)
1897 * NCFIV1: Nb of Coeff. of curves stored in CRBIV1. */
1898 /* CRBIV1: Table of coeffs of the approximation of iso-V0 and its */
1899 /* derivatives till order IORDRV. */
1900 /* NCFIV2: Nb of Coeff. of curves stored in CRBIV2. */
1901 /* CRBIV2: Table of coeffs of approximation of iso-V1 and its */
1902 /* derivatives till order IORDRV. */
1903 /* VHERMT: Coeff. of Hermit polynoms of order IORDRV. */
1904 /* PATJAC: Table of coefficients of the polynom P(u,v) of approximation */
1905 /* of F(u,v) WITHOUT taking into account the constraints. */
1907 /* OUTPUT ARGUMENTS : */
1908 /* ------------------- */
1909 /* PATJAC: Table of coefficients of the polynom P(u,v) by approximation */
1910 /* of F(u,v) WITH taking into account of constraints. */
1915 /* **********************************************************************
1917 /* Name of the routine */
1919 /* --------------------------- Initialisations --------------------------
1922 /* Parameter adjustments */
1923 patjac_dim1 = *mxujac + 1;
1924 patjac_dim2 = *mxvjac + 1;
1925 patjac_offset = patjac_dim1 * patjac_dim2;
1926 patjac -= patjac_offset;
1927 vhermt_dim1 = (*iordrv << 1) + 2;
1928 vhermt_offset = vhermt_dim1;
1929 vhermt -= vhermt_offset;
1932 crbiv2_dim1 = *nclimu;
1933 crbiv2_dim2 = *ndimen;
1934 crbiv2_offset = crbiv2_dim1 * (crbiv2_dim2 + 1);
1935 crbiv2 -= crbiv2_offset;
1936 crbiv1_dim1 = *nclimu;
1937 crbiv1_dim2 = *ndimen;
1938 crbiv1_offset = crbiv1_dim1 * (crbiv1_dim2 + 1);
1939 crbiv1 -= crbiv1_offset;
1942 ldbg = AdvApp2Var_SysBase::mnfndeb_() >= 3;
1944 AdvApp2Var_SysBase::mgenmsg_("MMA2AC2", 7L);
1947 /* ------------ ADDING of coeff by u of curves, by v of Hermit --------
1951 for (ii = 1; ii <= i__1; ++ii) {
1952 ndgv1 = ncfiv1[ii] - 1;
1953 ndgv2 = ncfiv2[ii] - 1;
1955 for (nd = 1; nd <= i__2; ++nd) {
1956 i__3 = (*iordrv << 1) + 1;
1957 for (jj = 0; jj <= i__3; ++jj) {
1958 bid1 = vhermt[jj + ((ii << 1) - 1) * vhermt_dim1];
1960 for (kk = 0; kk <= i__4; ++kk) {
1961 patjac[kk + (jj + nd * patjac_dim2) * patjac_dim1] +=
1962 bid1 * crbiv1[kk + (nd + ii * crbiv1_dim2) *
1966 bid2 = vhermt[jj + (ii << 1) * vhermt_dim1];
1968 for (kk = 0; kk <= i__4; ++kk) {
1969 patjac[kk + (jj + nd * patjac_dim2) * patjac_dim1] +=
1970 bid2 * crbiv2[kk + (nd + ii * crbiv2_dim2) *
1981 /* ------------------------------ The end -------------------------------
1985 AdvApp2Var_SysBase::mgsomsg_("MMA2AC2", 7L);
1991 //=======================================================================
1992 //function : mma2ac3_
1994 //=======================================================================
1995 int AdvApp2Var_ApproxF2var::mma2ac3_(const integer *ndimen,
1996 const integer *mxujac,
1997 const integer *mxvjac,
1998 const integer *iordru,
1999 const integer *nclimv,
2000 const integer *ncfiu1,
2001 const doublereal * crbiu1,
2002 const integer *ncfiu2,
2003 const doublereal *crbiu2,
2004 const doublereal *uhermt,
2008 /* System generated locals */
2009 integer crbiu1_dim1, crbiu1_dim2, crbiu1_offset, crbiu2_dim1, crbiu2_dim2,
2010 crbiu2_offset, patjac_dim1, patjac_dim2, patjac_offset,
2011 uhermt_dim1, uhermt_offset, i__1, i__2, i__3, i__4;
2013 /* Local variables */
2015 integer ndgu1, ndgu2, ii, jj, nd, kk;
2016 doublereal bid1, bid2;
2018 /* **********************************************************************
2023 /* Ajout des polynomes de contraintes */
2027 /* FONCTION,APPROXIMATION,COEFFICIENT,POLYNOME */
2029 /* INPUT ARGUMENTS : */
2030 /* ------------------ */
2031 /* NDIMEN: Dimension of the space. */
2032 /* MXUJAC: Max degree of the polynom of approximation by U. The */
2033 /* representation in the orthogonal base starts from degree */
2034 /* 0 to degree MXUJAC-2*(IORDRU+1). The polynomial base is the */
2035 /* base of Jacobi of order -1 (Legendre), 0, 1 or 2 */
2036 /* MXVJAC: Max degree of the polynom of approximation by V. The */
2037 /* representation in the orthogonal base starts from degree */
2038 /* 0 to degree MXUJAC-2*(IORDRU+1). The polynomial base is the */
2039 /* base of Jacobi of order -1 (Legendre), 0, 1 or 2 */
2040 /* IORDRU: Order of the base of Jacobi (-1,0,1 or 2) by U. Corresponds */
2041 /* to the step of constraints: C0, C1 or C2. */
2042 /* NCLIMV LIMIT nb of coeff by v of the solution P(u,v)
2043 * NCFIU1: Nb of Coeff. of curves stored in CRBIU1. */
2044 /* CRBIU1: Table of coeffs of the approximation of iso-U0 and its */
2045 /* derivatives till order IORDRU. */
2046 /* NCFIU2: Nb of Coeff. of curves stored in CRBIU2. */
2047 /* CRBIU2: Table of coeffs of approximation of iso-U1 and its */
2048 /* derivatives till order IORDRU */
2049 /* UHERMT: Coeff. of Hermit polynoms of order IORDRU. */
2050 /* PATJAC: Table of coefficients of the polynom P(u,v) of approximation */
2051 /* of F(u,v) WITHOUT taking into account the constraints. */
2053 /* OUTPUT ARGUMENTS : */
2054 /* ------------------- */
2055 /* PATJAC: Table of coefficients of the polynom P(u,v) by approximation */
2056 /* of F(u,v) WITH taking into account of constraints. */
2060 /* **********************************************************************
2062 /* The name of the routine */
2064 /* --------------------------- Initializations --------------------------
2067 /* Parameter adjustments */
2068 patjac_dim1 = *mxujac + 1;
2069 patjac_dim2 = *mxvjac + 1;
2070 patjac_offset = patjac_dim1 * patjac_dim2;
2071 patjac -= patjac_offset;
2072 uhermt_dim1 = (*iordru << 1) + 2;
2073 uhermt_offset = uhermt_dim1;
2074 uhermt -= uhermt_offset;
2077 crbiu2_dim1 = *nclimv;
2078 crbiu2_dim2 = *ndimen;
2079 crbiu2_offset = crbiu2_dim1 * (crbiu2_dim2 + 1);
2080 crbiu2 -= crbiu2_offset;
2081 crbiu1_dim1 = *nclimv;
2082 crbiu1_dim2 = *ndimen;
2083 crbiu1_offset = crbiu1_dim1 * (crbiu1_dim2 + 1);
2084 crbiu1 -= crbiu1_offset;
2087 ldbg = AdvApp2Var_SysBase::mnfndeb_() >= 3;
2089 AdvApp2Var_SysBase::mgenmsg_("MMA2AC3", 7L);
2092 /* ------------ ADDING of coeff by u of curves, by v of Hermit --------
2096 for (ii = 1; ii <= i__1; ++ii) {
2097 ndgu1 = ncfiu1[ii] - 1;
2098 ndgu2 = ncfiu2[ii] - 1;
2100 for (nd = 1; nd <= i__2; ++nd) {
2102 for (jj = 0; jj <= i__3; ++jj) {
2103 bid1 = crbiu1[jj + (nd + ii * crbiu1_dim2) * crbiu1_dim1];
2104 i__4 = (*iordru << 1) + 1;
2105 for (kk = 0; kk <= i__4; ++kk) {
2106 patjac[kk + (jj + nd * patjac_dim2) * patjac_dim1] +=
2107 bid1 * uhermt[kk + ((ii << 1) - 1) * uhermt_dim1];
2113 for (jj = 0; jj <= i__3; ++jj) {
2114 bid2 = crbiu2[jj + (nd + ii * crbiu2_dim2) * crbiu2_dim1];
2115 i__4 = (*iordru << 1) + 1;
2116 for (kk = 0; kk <= i__4; ++kk) {
2117 patjac[kk + (jj + nd * patjac_dim2) * patjac_dim1] +=
2118 bid2 * uhermt[kk + (ii << 1) * uhermt_dim1];
2129 /* ------------------------------ The end -------------------------------
2133 AdvApp2Var_SysBase::mgsomsg_("MMA2AC3", 7L);
2138 //=======================================================================
2139 //function : mma2can_
2141 //=======================================================================
2142 int AdvApp2Var_ApproxF2var::mma2can_(const integer *ncfmxu,
2143 const integer *ncfmxv,
2144 const integer *ndimen,
2145 const integer *iordru,
2146 const integer *iordrv,
2147 const integer *ncoefu,
2148 const integer *ncoefv,
2149 const doublereal *patjac,
2155 /* System generated locals */
2156 integer patjac_dim1, patjac_dim2, patjac_offset, patcan_dim1, patcan_dim2,
2157 patcan_offset, i__1, i__2;
2159 /* Local variables */
2161 integer ilon1, ilon2, ii, nd;
2163 /* **********************************************************************
2168 /* Change of Jacobi base to canonical (-1,1) and writing in a greater */
2173 /* ALL,AB_SPECIFI,CARREAU&,CONVERSION,JACOBI,CANNONIQUE,&CARREAU */
2175 /* INPUT ARGUMENTS : */
2176 /* -------------------- */
2177 /* NCFMXU: Dimension by U of resulting table PATCAN */
2178 /* NCFMXV: Dimension by V of resulting table PATCAN */
2179 /* NDIMEN: Dimension of the workspace. */
2180 /* IORDRU: Order of constraint by U */
2181 /* IORDRV: Order of constraint by V. */
2182 /* NCOEFU: Nb of coeff by U of square PATJAC */
2183 /* NCOEFV: Nb of coeff by V of square PATJAC */
2184 /* PATJAC: Square in the base of Jacobi of order IORDRU by U and */
2187 /* OUTPUT ARGUMENTS : */
2188 /* --------------------- */
2189 /* PATAUX: Auxiliary Table. */
2190 /* PATCAN: Table of coefficients in the canonic base. */
2191 /* IERCOD: Error code. */
2192 /* = 0, everything goes well, and all things are equal. */
2193 /* = 1, the program refuses to process with incorrect input arguments */
2196 /* COMMONS USED : */
2197 /* ------------------ */
2199 /* REFERENCES CALLED : */
2200 /* --------------------- */
2202 /* DESCRIPTION/NOTES/LIMITATIONS : */
2203 /* ----------------------------------- */
2205 /* **********************************************************************
2209 /* Parameter adjustments */
2210 patcan_dim1 = *ncfmxu;
2211 patcan_dim2 = *ncfmxv;
2212 patcan_offset = patcan_dim1 * (patcan_dim2 + 1) + 1;
2213 patcan -= patcan_offset;
2215 patjac_dim1 = *ncoefu;
2216 patjac_dim2 = *ncoefv;
2217 patjac_offset = patjac_dim1 * (patjac_dim2 + 1) + 1;
2218 patjac -= patjac_offset;
2221 ldbg = AdvApp2Var_SysBase::mnfndeb_() >= 2;
2223 AdvApp2Var_SysBase::mgenmsg_("MMA2CAN", 7L);
2227 if (*iordru < -1 || *iordru > 2) {
2230 if (*iordrv < -1 || *iordrv > 2) {
2233 if (*ncoefu > *ncfmxu || *ncoefv > *ncfmxv) {
2237 /* --> Pass to canonic base (-1,1) */
2238 mmjacpt_(ndimen, ncoefu, ncoefv, iordru, iordrv, &patjac[patjac_offset], &
2239 pataux[1], &patcan[patcan_offset]);
2241 /* --> Write all in a greater table */
2242 AdvApp2Var_MathBase::mmfmca8_(ncoefu,
2248 &patcan[patcan_offset],
2249 &patcan[patcan_offset]);
2251 /* --> Complete with zeros the resulting table. */
2252 ilon1 = *ncfmxu - *ncoefu;
2253 ilon2 = *ncfmxu * (*ncfmxv - *ncoefv);
2255 for (nd = 1; nd <= i__1; ++nd) {
2258 for (ii = 1; ii <= i__2; ++ii) {
2259 AdvApp2Var_SysBase::mvriraz_(&ilon1,
2260 &patcan[*ncoefu + 1 + (ii + nd * patcan_dim2) * patcan_dim1]);
2265 AdvApp2Var_SysBase::mvriraz_(&ilon2,
2266 &patcan[(*ncoefv + 1 + nd * patcan_dim2) * patcan_dim1 + 1]);
2273 /* ----------------------
2281 AdvApp2Var_SysBase::maermsg_("MMA2CAN", iercod, 7L);
2283 AdvApp2Var_SysBase::mgsomsg_("MMA2CAN", 7L);
2288 //=======================================================================
2289 //function : mma2cd1_
2291 //=======================================================================
2292 int mma2cd1_(integer *ndimen,
2315 /* System generated locals */
2316 integer contr1_dim1, contr1_dim2, contr1_offset, contr2_dim1, contr2_dim2,
2317 contr2_offset, contr3_dim1, contr3_dim2, contr3_offset,
2318 contr4_dim1, contr4_dim2, contr4_offset, uhermt_dim1,
2319 uhermt_offset, vhermt_dim1, vhermt_offset, fpntbu_dim1,
2320 fpntbu_offset, fpntbv_dim1, fpntbv_offset, sosotb_dim1,
2321 sosotb_dim2, sosotb_offset, diditb_dim1, diditb_dim2,
2322 diditb_offset, soditb_dim1, soditb_dim2, soditb_offset,
2323 disotb_dim1, disotb_dim2, disotb_offset, i__1, i__2, i__3, i__4,
2326 /* Local variables */
2327 integer ncfhu, ncfhv, nuroo, nvroo, nd, ii, jj, kk, ll, ibb, kkm,
2329 doublereal bid1, bid2, bid3, bid4;
2330 doublereal diu1, diu2, div1, div2, sou1, sou2, sov1, sov2;
2332 /* **********************************************************************
2337 /* Discretisation on the parameters of polynoms of interpolation */
2338 /* of constraints at the corners of order IORDRE. */
2342 /* TOUS, AB_SPECIFI::CONTRAINTE&, DISCRETISATION, &POINT */
2344 /* INPUT ARGUMENTS : */
2345 /* ------------------ */
2346 /* NDIMEN: Dimension of the space. */
2347 /* NBPNTU: Nb of INTERNAL parameters of discretisation by U. */
2348 /* This is also the nb of root of Legendre polynom where discretization is done. */
2349 /* UROOTL: Table of parameters of discretisation ON (-1,1) by U.
2351 /* NBPNTV: Nb of INTERNAL parameters of discretisation by V. */
2352 /* This is also the nb of root of Legendre polynom where discretization is done. */
2353 /* VROOTL: Table of discretization parameters on (-1,1) by V.
2354 /* IORDRU: Order of constraint imposed at the extremities of iso-V */
2355 /* = 0, calculate the extremities of iso-V */
2356 /* = 1, calculate, additionally, the 1st derivative in the direction of iso-V */
2357 /* = 2, calculate, additionally, the 2nd derivative in the direction of iso-V */
2358 /* IORDRV: Order of constraint imposed at the extremities of iso-U */
2359 /* = 0, calculate the extremities of iso-U */
2360 /* = 1, calculate, additionally, the 1st derivative in the direction of iso-U */
2361 /* = 2, calculate, additionally, the 2nd derivative in the direction of iso-U */
2362 /* CONTR1: Contains, if IORDRU and IORDRV>=0, the values at the */
2363 /* extremities of F(U0,V0) and its derivatives. */
2364 /* CONTR2: Contains, if IORDRU and IORDRV>=0, the values at the */
2365 /* extremities of F(U1,V0) and its derivatives. */
2366 /* CONTR3: Contains, if IORDRU and IORDRV>=0, the values at the */
2367 /* extremities of F(U0,V1) and its derivatives. */
2368 /* CONTR4: Contains, if IORDRU and IORDRV>=0, the values at the */
2369 /* extremities of F(U1,V1) and its derivatives. */
2370 /* SOSOTB: Preinitialized table (input/output argument). */
2371 /* DISOTB: Preinitialized table (input/output argument). */
2372 /* SODITB: Preinitialized table (input/output argument). */
2373 /* DIDITB: Preinitialized table (input/output argument) */
2375 /* OUTPUT ARGUMENTS : */
2376 /* ------------------- */
2377 /* FPNTBU: Auxiliary table. */
2378 /* FPNTBV: Auxiliary table. */
2379 /* UHERMT: Table of 2*(IORDRU+1) coeff. of 2*(IORDRU+1) polynoms of Hermite. */
2380 /* VHERMT: Table of 2*(IORDRV+1) coeff. of 2*(IORDRV+1) polynoms of Hermite. */
2381 /* SOSOTB: Table where the terms of constraints are added */
2382 /* C(ui,vj) + C(ui,-vj) + C(-ui,vj) + C(-ui,-vj) */
2383 /* with ui and vj positive roots of the Legendre polynom */
2384 /* of degree NBPNTU and NBPNTV respectively. */
2385 /* DISOTB: Table where the terms of constraints are added */
2386 /* C(ui,vj) + C(ui,-vj) - C(-ui,vj) - C(-ui,-vj) */
2387 /* with ui and vj positive roots of the polynom of Legendre */
2388 /* of degree NBPNTU and NBPNTV respectively. */
2389 /* SODITB: Table where the terms of constraints are added */
2390 /* C(ui,vj) - C(ui,-vj) + C(-ui,vj) - C(-ui,-vj) */
2391 /* with ui and vj positive roots of the polynom of Legendre */
2392 /* of degree NBPNTU and NBPNTV respectively. */
2393 /* DIDITB: Table where the terms of constraints are added */
2394 /* C(ui,vj) - C(ui,-vj) - C(-ui,vj) + C(-ui,-vj) */
2395 /* with ui and vj positive roots of the polynom of Legendre */
2396 /* of degree NBPNTU and NBPNTV respectively. */
2398 /* COMMONS USED : */
2399 /* ---------------- */
2401 /* REFERENCES CALLED : */
2402 /* ----------------------- */
2404 /* DESCRIPTION/NOTES/LIMITATIONS : */
2405 /* ----------------------------------- */
2408 /* **********************************************************************
2411 /* Name of the routine */
2414 /* Parameter adjustments */
2416 diditb_dim1 = *nbpntu / 2 + 1;
2417 diditb_dim2 = *nbpntv / 2 + 1;
2418 diditb_offset = diditb_dim1 * diditb_dim2;
2419 diditb -= diditb_offset;
2420 disotb_dim1 = *nbpntu / 2;
2421 disotb_dim2 = *nbpntv / 2;
2422 disotb_offset = disotb_dim1 * (disotb_dim2 + 1) + 1;
2423 disotb -= disotb_offset;
2424 soditb_dim1 = *nbpntu / 2;
2425 soditb_dim2 = *nbpntv / 2;
2426 soditb_offset = soditb_dim1 * (soditb_dim2 + 1) + 1;
2427 soditb -= soditb_offset;
2428 sosotb_dim1 = *nbpntu / 2 + 1;
2429 sosotb_dim2 = *nbpntv / 2 + 1;
2430 sosotb_offset = sosotb_dim1 * sosotb_dim2;
2431 sosotb -= sosotb_offset;
2433 uhermt_dim1 = (*iordru << 1) + 2;
2434 uhermt_offset = uhermt_dim1;
2435 uhermt -= uhermt_offset;
2436 fpntbu_dim1 = *nbpntu;
2437 fpntbu_offset = fpntbu_dim1 + 1;
2438 fpntbu -= fpntbu_offset;
2439 vhermt_dim1 = (*iordrv << 1) + 2;
2440 vhermt_offset = vhermt_dim1;
2441 vhermt -= vhermt_offset;
2442 fpntbv_dim1 = *nbpntv;
2443 fpntbv_offset = fpntbv_dim1 + 1;
2444 fpntbv -= fpntbv_offset;
2445 contr4_dim1 = *ndimen;
2446 contr4_dim2 = *iordru + 2;
2447 contr4_offset = contr4_dim1 * (contr4_dim2 + 1) + 1;
2448 contr4 -= contr4_offset;
2449 contr3_dim1 = *ndimen;
2450 contr3_dim2 = *iordru + 2;
2451 contr3_offset = contr3_dim1 * (contr3_dim2 + 1) + 1;
2452 contr3 -= contr3_offset;
2453 contr2_dim1 = *ndimen;
2454 contr2_dim2 = *iordru + 2;
2455 contr2_offset = contr2_dim1 * (contr2_dim2 + 1) + 1;
2456 contr2 -= contr2_offset;
2457 contr1_dim1 = *ndimen;
2458 contr1_dim2 = *iordru + 2;
2459 contr1_offset = contr1_dim1 * (contr1_dim2 + 1) + 1;
2460 contr1 -= contr1_offset;
2463 ibb = AdvApp2Var_SysBase::mnfndeb_();
2465 AdvApp2Var_SysBase::mgenmsg_("MMA2CD1", 7L);
2468 /* ------------------- Discretisation of Hermite polynoms -----------
2471 ncfhu = (*iordru + 1) << 1;
2473 for (ii = 1; ii <= i__1; ++ii) {
2475 for (ll = 1; ll <= i__2; ++ll) {
2476 AdvApp2Var_MathBase::mmmpocur_(&ncfhu, &c__1, &ncfhu, &uhermt[ii * uhermt_dim1], &
2477 urootl[ll], &fpntbu[ll + ii * fpntbu_dim1]);
2482 ncfhv = (*iordrv + 1) << 1;
2484 for (jj = 1; jj <= i__1; ++jj) {
2486 for (kk = 1; kk <= i__2; ++kk) {
2487 AdvApp2Var_MathBase::mmmpocur_(&ncfhv, &c__1, &ncfhv, &vhermt[jj * vhermt_dim1], &
2488 vrootl[kk], &fpntbv[kk + jj * fpntbv_dim1]);
2494 /* ---- The discretizations of polynoms of constraints are subtracted ----
2497 nuroo = *nbpntu / 2;
2498 nvroo = *nbpntv / 2;
2500 for (nd = 1; nd <= i__1; ++nd) {
2503 for (jj = 1; jj <= i__2; ++jj) {
2505 for (ii = 1; ii <= i__3; ++ii) {
2506 bid1 = contr1[nd + (ii + jj * contr1_dim2) * contr1_dim1];
2507 bid2 = contr2[nd + (ii + jj * contr2_dim2) * contr2_dim1];
2508 bid3 = contr3[nd + (ii + jj * contr3_dim2) * contr3_dim1];
2509 bid4 = contr4[nd + (ii + jj * contr4_dim2) * contr4_dim1];
2512 for (kk = 1; kk <= i__4; ++kk) {
2513 kkp = (*nbpntv + 1) / 2 + kk;
2514 kkm = nvroo - kk + 1;
2515 sov1 = fpntbv[kkp + ((jj << 1) - 1) * fpntbv_dim1] +
2516 fpntbv[kkm + ((jj << 1) - 1) * fpntbv_dim1];
2517 div1 = fpntbv[kkp + ((jj << 1) - 1) * fpntbv_dim1] -
2518 fpntbv[kkm + ((jj << 1) - 1) * fpntbv_dim1];
2519 sov2 = fpntbv[kkp + (jj << 1) * fpntbv_dim1] + fpntbv[kkm
2520 + (jj << 1) * fpntbv_dim1];
2521 div2 = fpntbv[kkp + (jj << 1) * fpntbv_dim1] - fpntbv[kkm
2522 + (jj << 1) * fpntbv_dim1];
2524 for (ll = 1; ll <= i__5; ++ll) {
2525 llp = (*nbpntu + 1) / 2 + ll;
2526 llm = nuroo - ll + 1;
2527 sou1 = fpntbu[llp + ((ii << 1) - 1) * fpntbu_dim1] +
2528 fpntbu[llm + ((ii << 1) - 1) * fpntbu_dim1];
2529 diu1 = fpntbu[llp + ((ii << 1) - 1) * fpntbu_dim1] -
2530 fpntbu[llm + ((ii << 1) - 1) * fpntbu_dim1];
2531 sou2 = fpntbu[llp + (ii << 1) * fpntbu_dim1] + fpntbu[
2532 llm + (ii << 1) * fpntbu_dim1];
2533 diu2 = fpntbu[llp + (ii << 1) * fpntbu_dim1] - fpntbu[
2534 llm + (ii << 1) * fpntbu_dim1];
2535 sosotb[ll + (kk + nd * sosotb_dim2) * sosotb_dim1] =
2536 sosotb[ll + (kk + nd * sosotb_dim2) *
2537 sosotb_dim1] - bid1 * sou1 * sov1 - bid2 *
2538 sou2 * sov1 - bid3 * sou1 * sov2 - bid4 *
2540 soditb[ll + (kk + nd * soditb_dim2) * soditb_dim1] =
2541 soditb[ll + (kk + nd * soditb_dim2) *
2542 soditb_dim1] - bid1 * sou1 * div1 - bid2 *
2543 sou2 * div1 - bid3 * sou1 * div2 - bid4 *
2545 disotb[ll + (kk + nd * disotb_dim2) * disotb_dim1] =
2546 disotb[ll + (kk + nd * disotb_dim2) *
2547 disotb_dim1] - bid1 * diu1 * sov1 - bid2 *
2548 diu2 * sov1 - bid3 * diu1 * sov2 - bid4 *
2550 diditb[ll + (kk + nd * diditb_dim2) * diditb_dim1] =
2551 diditb[ll + (kk + nd * diditb_dim2) *
2552 diditb_dim1] - bid1 * diu1 * div1 - bid2 *
2553 diu2 * div1 - bid3 * diu1 * div2 - bid4 *
2560 /* ------------ Case when the discretization is done only on the roots
2562 /* ---------- of Legendre polynom of uneven degree, 0 is root
2565 if (*nbpntu % 2 == 1) {
2566 sou1 = fpntbu[nuroo + 1 + ((ii << 1) - 1) * fpntbu_dim1];
2567 sou2 = fpntbu[nuroo + 1 + (ii << 1) * fpntbu_dim1];
2569 for (kk = 1; kk <= i__4; ++kk) {
2570 kkp = (*nbpntv + 1) / 2 + kk;
2571 kkm = nvroo - kk + 1;
2572 sov1 = fpntbv[kkp + ((jj << 1) - 1) * fpntbv_dim1] +
2573 fpntbv[kkm + ((jj << 1) - 1) * fpntbv_dim1];
2574 div1 = fpntbv[kkp + ((jj << 1) - 1) * fpntbv_dim1] -
2575 fpntbv[kkm + ((jj << 1) - 1) * fpntbv_dim1];
2576 sov2 = fpntbv[kkp + (jj << 1) * fpntbv_dim1] + fpntbv[
2577 kkm + (jj << 1) * fpntbv_dim1];
2578 div2 = fpntbv[kkp + (jj << 1) * fpntbv_dim1] - fpntbv[
2579 kkm + (jj << 1) * fpntbv_dim1];
2580 sosotb[(kk + nd * sosotb_dim2) * sosotb_dim1] =
2581 sosotb[(kk + nd * sosotb_dim2) * sosotb_dim1]
2582 - bid1 * sou1 * sov1 - bid2 * sou2 * sov1 -
2583 bid3 * sou1 * sov2 - bid4 * sou2 * sov2;
2584 diditb[(kk + nd * diditb_dim2) * diditb_dim1] =
2585 diditb[(kk + nd * diditb_dim2) * diditb_dim1]
2586 - bid1 * sou1 * div1 - bid2 * sou2 * div1 -
2587 bid3 * sou1 * div2 - bid4 * sou2 * div2;
2592 if (*nbpntv % 2 == 1) {
2593 sov1 = fpntbv[nvroo + 1 + ((jj << 1) - 1) * fpntbv_dim1];
2594 sov2 = fpntbv[nvroo + 1 + (jj << 1) * fpntbv_dim1];
2596 for (ll = 1; ll <= i__4; ++ll) {
2597 llp = (*nbpntu + 1) / 2 + ll;
2598 llm = nuroo - ll + 1;
2599 sou1 = fpntbu[llp + ((ii << 1) - 1) * fpntbu_dim1] +
2600 fpntbu[llm + ((ii << 1) - 1) * fpntbu_dim1];
2601 diu1 = fpntbu[llp + ((ii << 1) - 1) * fpntbu_dim1] -
2602 fpntbu[llm + ((ii << 1) - 1) * fpntbu_dim1];
2603 sou2 = fpntbu[llp + (ii << 1) * fpntbu_dim1] + fpntbu[
2604 llm + (ii << 1) * fpntbu_dim1];
2605 diu2 = fpntbu[llp + (ii << 1) * fpntbu_dim1] - fpntbu[
2606 llm + (ii << 1) * fpntbu_dim1];
2607 sosotb[ll + nd * sosotb_dim2 * sosotb_dim1] = sosotb[
2608 ll + nd * sosotb_dim2 * sosotb_dim1] - bid1 *
2609 sou1 * sov1 - bid2 * sou2 * sov1 - bid3 *
2610 sou1 * sov2 - bid4 * sou2 * sov2;
2611 diditb[ll + nd * diditb_dim2 * diditb_dim1] = diditb[
2612 ll + nd * diditb_dim2 * diditb_dim1] - bid1 *
2613 diu1 * sov1 - bid2 * diu2 * sov1 - bid3 *
2614 diu1 * sov2 - bid4 * diu2 * sov2;
2619 if (*nbpntu % 2 == 1 && *nbpntv % 2 == 1) {
2620 sou1 = fpntbu[nuroo + 1 + ((ii << 1) - 1) * fpntbu_dim1];
2621 sou2 = fpntbu[nuroo + 1 + (ii << 1) * fpntbu_dim1];
2622 sov1 = fpntbv[nvroo + 1 + ((jj << 1) - 1) * fpntbv_dim1];
2623 sov2 = fpntbv[nvroo + 1 + (jj << 1) * fpntbv_dim1];
2624 sosotb[nd * sosotb_dim2 * sosotb_dim1] = sosotb[nd *
2625 sosotb_dim2 * sosotb_dim1] - bid1 * sou1 * sov1 -
2626 bid2 * sou2 * sov1 - bid3 * sou1 * sov2 - bid4 *
2628 diditb[nd * diditb_dim2 * diditb_dim1] = diditb[nd *
2629 diditb_dim2 * diditb_dim1] - bid1 * sou1 * sov1 -
2630 bid2 * sou2 * sov1 - bid3 * sou1 * sov2 - bid4 *
2642 /* ------------------------------ The End -------------------------------
2647 AdvApp2Var_SysBase::mgsomsg_("MMA2CD1", 7L);
2652 //=======================================================================
2653 //function : mma2cd2_
2655 //=======================================================================
2656 int mma2cd2_(integer *ndimen,
2674 /* System generated locals */
2675 integer sotbv1_dim1, sotbv1_dim2, sotbv1_offset, sotbv2_dim1, sotbv2_dim2,
2676 sotbv2_offset, ditbv1_dim1, ditbv1_dim2, ditbv1_offset,
2677 ditbv2_dim1, ditbv2_dim2, ditbv2_offset, fpntab_dim1,
2678 fpntab_offset, vhermt_dim1, vhermt_offset, sosotb_dim1,
2679 sosotb_dim2, sosotb_offset, diditb_dim1, diditb_dim2,
2680 diditb_offset, soditb_dim1, soditb_dim2, soditb_offset,
2681 disotb_dim1, disotb_dim2, disotb_offset, i__1, i__2, i__3, i__4;
2683 /* Local variables */
2684 integer ncfhv, nuroo, nvroo, ii, nd, jj, kk, ibb, jjm, jjp;
2685 doublereal bid1, bid2, bid3, bid4;
2687 /* **********************************************************************
2691 /* Discretisation on the parameters of polynoms of interpolation */
2692 /* of constraints on 2 borders iso-V of order IORDRV. */
2697 /* TOUS, AB_SPECIFI::CONTRAINTE&, DISCRETISATION, &POINT */
2701 /* INPUT ARGUMENTS : */
2702 /* ------------------ */
2703 /* NDIMEN: Dimension of the space. */
2704 /* NBPNTU: Nb of INTERNAL parameters of discretisation by U. */
2705 /* This is also the nb of root of Legendre polynom where discretization is done. */
2706 /* UROOTL: Table of parameters of discretisation ON (-1,1) by U.
2708 /* NBPNTV: Nb of INTERNAL parameters of discretisation by V. */
2709 /* This is also the nb of root of Legendre polynom where discretization is done. */
2710 /* VROOTL: Table of discretization parameters on (-1,1) by V.
2711 /* IORDRV: Order of constraint imposed at the extremities of iso-V */
2712 /* = 0, calculate the extremities of iso-V */
2713 /* = 1, calculate, additionally, the 1st derivative in the direction of iso-V */
2714 /* = 2, calculate, additionally, the 2nd derivative in the direction of iso-V */
2715 /* SOTBV1: Table of NBPNTV/2 sums of 2 index points */
2716 /* NBPNTV-II+1 and II, for II = 1, NBPNTV/2 on iso-V0. */
2717 /* SOTBV2: Table of NBPNTV/2 sums of 2 index points */
2718 /* NBPNTV-II+1 and II, for II = 1, NBPNTV/2 on iso-V1. */
2719 /* DITBV1: Table of NBPNTV/2 differences of 2 index points */
2720 /* NBPNTV-II+1 and II, for II = 1, NBPNTV/2 on iso-V0. */
2721 /* DITBV2: Table of NBPNTV/2 differences of 2 index points */
2722 /* NBPNTV-II+1 and II, for II = 1, NBPNTV/2 on iso-V1. */
2723 /* SOSOTB: Preinitialized table (input/output argument). */
2724 /* DISOTB: Preinitialized table (input/output argument). */
2725 /* SODITB: Preinitialized table (input/output argument). */
2726 /* DIDITB: Preinitialized table (input/output argument) */
2728 /* OUTPUT ARGUMENTS : */
2729 /* ------------------- */
2730 /* FPNTAB: Auxiliary table. */
2731 /* VHERMT: Table of 2*(IORDRV+1) coeff. of 2*(IORDRV+1) polynoms of Hermite. */
2732 /* SOSOTB: Table where the terms of constraints are added */
2733 /* C(ui,vj) + C(ui,-vj) + C(-ui,vj) + C(-ui,-vj) */
2734 /* with ui and vj positive roots of the Legendre polynom */
2735 /* of degree NBPNTU and NBPNTV respectively. */
2736 /* DISOTB: Table where the terms of constraints are added */
2737 /* C(ui,vj) + C(ui,-vj) - C(-ui,vj) - C(-ui,-vj) */
2738 /* with ui and vj positive roots of the polynom of Legendre */
2739 /* of degree NBPNTU and NBPNTV respectively. */
2740 /* SODITB: Table where the terms of constraints are added */
2741 /* C(ui,vj) - C(ui,-vj) + C(-ui,vj) - C(-ui,-vj) */
2742 /* with ui and vj positive roots of the polynom of Legendre */
2743 /* of degree NBPNTU and NBPNTV respectively. */
2744 /* DIDITB: Table where the terms of constraints are added */
2745 /* C(ui,vj) - C(ui,-vj) - C(-ui,vj) + C(-ui,-vj) */
2746 /* with ui and vj positive roots of the polynom of Legendre */
2747 /* of degree NBPNTU and NBPNTV respectively. */
2749 /* COMMONS USED : */
2750 /* ---------------- */
2752 /* REFERENCES CALLED : */
2753 /* ----------------------- */
2755 /* DESCRIPTION/NOTES/LIMITATIONS : */
2756 /* ----------------------------------- */
2760 /* **********************************************************************
2763 /* Name of the routine */
2766 /* Parameter adjustments */
2767 diditb_dim1 = *nbpntu / 2 + 1;
2768 diditb_dim2 = *nbpntv / 2 + 1;
2769 diditb_offset = diditb_dim1 * diditb_dim2;
2770 diditb -= diditb_offset;
2771 disotb_dim1 = *nbpntu / 2;
2772 disotb_dim2 = *nbpntv / 2;
2773 disotb_offset = disotb_dim1 * (disotb_dim2 + 1) + 1;
2774 disotb -= disotb_offset;
2775 soditb_dim1 = *nbpntu / 2;
2776 soditb_dim2 = *nbpntv / 2;
2777 soditb_offset = soditb_dim1 * (soditb_dim2 + 1) + 1;
2778 soditb -= soditb_offset;
2779 sosotb_dim1 = *nbpntu / 2 + 1;
2780 sosotb_dim2 = *nbpntv / 2 + 1;
2781 sosotb_offset = sosotb_dim1 * sosotb_dim2;
2782 sosotb -= sosotb_offset;
2784 vhermt_dim1 = (*iordrv << 1) + 2;
2785 vhermt_offset = vhermt_dim1;
2786 vhermt -= vhermt_offset;
2787 fpntab_dim1 = *nbpntv;
2788 fpntab_offset = fpntab_dim1 + 1;
2789 fpntab -= fpntab_offset;
2790 ditbv2_dim1 = *nbpntu / 2 + 1;
2791 ditbv2_dim2 = *ndimen;
2792 ditbv2_offset = ditbv2_dim1 * (ditbv2_dim2 + 1);
2793 ditbv2 -= ditbv2_offset;
2794 ditbv1_dim1 = *nbpntu / 2 + 1;
2795 ditbv1_dim2 = *ndimen;
2796 ditbv1_offset = ditbv1_dim1 * (ditbv1_dim2 + 1);
2797 ditbv1 -= ditbv1_offset;
2798 sotbv2_dim1 = *nbpntu / 2 + 1;
2799 sotbv2_dim2 = *ndimen;
2800 sotbv2_offset = sotbv2_dim1 * (sotbv2_dim2 + 1);
2801 sotbv2 -= sotbv2_offset;
2802 sotbv1_dim1 = *nbpntu / 2 + 1;
2803 sotbv1_dim2 = *ndimen;
2804 sotbv1_offset = sotbv1_dim1 * (sotbv1_dim2 + 1);
2805 sotbv1 -= sotbv1_offset;
2808 ibb = AdvApp2Var_SysBase::mnfndeb_();
2810 AdvApp2Var_SysBase::mgenmsg_("MMA2CD2", 7L);
2813 /* ------------------- Discretization of Hermit polynoms -----------
2816 ncfhv = (*iordrv + 1) << 1;
2818 for (ii = 1; ii <= i__1; ++ii) {
2820 for (jj = 1; jj <= i__2; ++jj) {
2821 AdvApp2Var_MathBase::mmmpocur_(&ncfhv, &c__1, &ncfhv, &vhermt[ii * vhermt_dim1], &
2822 vrootl[jj], &fpntab[jj + ii * fpntab_dim1]);
2828 /* ---- The discretizations of polynoms of constraints are subtracted ----
2831 nuroo = *nbpntu / 2;
2832 nvroo = *nbpntv / 2;
2835 for (nd = 1; nd <= i__1; ++nd) {
2837 for (ii = 1; ii <= i__2; ++ii) {
2840 for (kk = 1; kk <= i__3; ++kk) {
2841 bid1 = sotbv1[kk + (nd + ii * sotbv1_dim2) * sotbv1_dim1];
2842 bid2 = sotbv2[kk + (nd + ii * sotbv2_dim2) * sotbv2_dim1];
2843 bid3 = ditbv1[kk + (nd + ii * ditbv1_dim2) * ditbv1_dim1];
2844 bid4 = ditbv2[kk + (nd + ii * ditbv2_dim2) * ditbv2_dim1];
2846 for (jj = 1; jj <= i__4; ++jj) {
2847 jjp = (*nbpntv + 1) / 2 + jj;
2848 jjm = nvroo - jj + 1;
2849 sosotb[kk + (jj + nd * sosotb_dim2) * sosotb_dim1] =
2850 sosotb[kk + (jj + nd * sosotb_dim2) * sosotb_dim1]
2851 - bid1 * (fpntab[jjp + ((ii << 1) - 1) *
2852 fpntab_dim1] + fpntab[jjm + ((ii << 1) - 1) *
2853 fpntab_dim1]) - bid2 * (fpntab[jjp + (ii << 1) *
2854 fpntab_dim1] + fpntab[jjm + (ii << 1) *
2856 disotb[kk + (jj + nd * disotb_dim2) * disotb_dim1] =
2857 disotb[kk + (jj + nd * disotb_dim2) * disotb_dim1]
2858 - bid3 * (fpntab[jjp + ((ii << 1) - 1) *
2859 fpntab_dim1] + fpntab[jjm + ((ii << 1) - 1) *
2860 fpntab_dim1]) - bid4 * (fpntab[jjp + (ii << 1) *
2861 fpntab_dim1] + fpntab[jjm + (ii << 1) *
2863 soditb[kk + (jj + nd * soditb_dim2) * soditb_dim1] =
2864 soditb[kk + (jj + nd * soditb_dim2) * soditb_dim1]
2865 - bid1 * (fpntab[jjp + ((ii << 1) - 1) *
2866 fpntab_dim1] - fpntab[jjm + ((ii << 1) - 1) *
2867 fpntab_dim1]) - bid2 * (fpntab[jjp + (ii << 1) *
2868 fpntab_dim1] - fpntab[jjm + (ii << 1) *
2870 diditb[kk + (jj + nd * diditb_dim2) * diditb_dim1] =
2871 diditb[kk + (jj + nd * diditb_dim2) * diditb_dim1]
2872 - bid3 * (fpntab[jjp + ((ii << 1) - 1) *
2873 fpntab_dim1] - fpntab[jjm + ((ii << 1) - 1) *
2874 fpntab_dim1]) - bid4 * (fpntab[jjp + (ii << 1) *
2875 fpntab_dim1] - fpntab[jjm + (ii << 1) *
2884 /* ------------ Case when the discretization is done only on the roots */
2885 /* ---------- of Legendre polynom of uneven degree, 0 is root */
2888 if (*nbpntv % 2 == 1) {
2890 for (ii = 1; ii <= i__2; ++ii) {
2892 for (kk = 1; kk <= i__3; ++kk) {
2893 bid1 = sotbv1[kk + (nd + ii * sotbv1_dim2) * sotbv1_dim1]
2894 * fpntab[nvroo + 1 + ((ii << 1) - 1) *
2895 fpntab_dim1] + sotbv2[kk + (nd + ii * sotbv2_dim2)
2896 * sotbv2_dim1] * fpntab[nvroo + 1 + (ii << 1) *
2898 sosotb[kk + nd * sosotb_dim2 * sosotb_dim1] -= bid1;
2899 bid2 = ditbv1[kk + (nd + ii * ditbv1_dim2) * ditbv1_dim1]
2900 * fpntab[nvroo + 1 + ((ii << 1) - 1) *
2901 fpntab_dim1] + ditbv2[kk + (nd + ii * ditbv2_dim2)
2902 * ditbv2_dim1] * fpntab[nvroo + 1 + (ii << 1) *
2904 diditb[kk + nd * diditb_dim2 * diditb_dim1] -= bid2;
2911 if (*nbpntu % 2 == 1) {
2913 for (ii = 1; ii <= i__2; ++ii) {
2915 for (jj = 1; jj <= i__3; ++jj) {
2916 jjp = (*nbpntv + 1) / 2 + jj;
2917 jjm = nvroo - jj + 1;
2918 bid1 = sotbv1[(nd + ii * sotbv1_dim2) * sotbv1_dim1] * (
2919 fpntab[jjp + ((ii << 1) - 1) * fpntab_dim1] +
2920 fpntab[jjm + ((ii << 1) - 1) * fpntab_dim1]) +
2921 sotbv2[(nd + ii * sotbv2_dim2) * sotbv2_dim1] * (
2922 fpntab[jjp + (ii << 1) * fpntab_dim1] + fpntab[
2923 jjm + (ii << 1) * fpntab_dim1]);
2924 sosotb[(jj + nd * sosotb_dim2) * sosotb_dim1] -= bid1;
2925 bid2 = sotbv1[(nd + ii * sotbv1_dim2) * sotbv1_dim1] * (
2926 fpntab[jjp + ((ii << 1) - 1) * fpntab_dim1] -
2927 fpntab[jjm + ((ii << 1) - 1) * fpntab_dim1]) +
2928 sotbv2[(nd + ii * sotbv2_dim2) * sotbv2_dim1] * (
2929 fpntab[jjp + (ii << 1) * fpntab_dim1] - fpntab[
2930 jjm + (ii << 1) * fpntab_dim1]);
2931 diditb[jj + nd * diditb_dim2 * diditb_dim1] -= bid2;
2938 if (*nbpntu % 2 == 1 && *nbpntv % 2 == 1) {
2940 for (ii = 1; ii <= i__2; ++ii) {
2941 bid1 = sotbv1[(nd + ii * sotbv1_dim2) * sotbv1_dim1] * fpntab[
2942 nvroo + 1 + ((ii << 1) - 1) * fpntab_dim1] + sotbv2[(
2943 nd + ii * sotbv2_dim2) * sotbv2_dim1] * fpntab[nvroo
2944 + 1 + (ii << 1) * fpntab_dim1];
2945 sosotb[nd * sosotb_dim2 * sosotb_dim1] -= bid1;
2954 /* ------------------------------ The End -------------------------------
2959 AdvApp2Var_SysBase::mgsomsg_("MMA2CD2", 7L);
2964 //=======================================================================
2965 //function : mma2cd3_
2967 //=======================================================================
2968 int mma2cd3_(integer *ndimen,
2987 /* System generated locals */
2988 integer sotbu1_dim1, sotbu1_dim2, sotbu1_offset, sotbu2_dim1, sotbu2_dim2,
2989 sotbu2_offset, ditbu1_dim1, ditbu1_dim2, ditbu1_offset,
2990 ditbu2_dim1, ditbu2_dim2, ditbu2_offset, fpntab_dim1,
2991 fpntab_offset, uhermt_dim1, uhermt_offset, sosotb_dim1,
2992 sosotb_dim2, sosotb_offset, diditb_dim1, diditb_dim2,
2993 diditb_offset, soditb_dim1, soditb_dim2, soditb_offset,
2994 disotb_dim1, disotb_dim2, disotb_offset, i__1, i__2, i__3, i__4;
2996 /* Local variables */
2997 integer ncfhu, nuroo, nvroo, ii, nd, jj, kk, ibb, kkm, kkp;
2998 doublereal bid1, bid2, bid3, bid4;
3000 /* **********************************************************************
3004 /* Discretisation on the parameters of polynoms of interpolation */
3005 /* of constraints on 2 borders iso-U of order IORDRU. */
3010 /* TOUS, AB_SPECIFI::CONTRAINTE&, DISCRETISATION, &POINT */
3012 /* INPUT ARGUMENTS : */
3013 /* ------------------ */
3014 /* NDIMEN: Dimension of the space. */
3015 /* NBPNTU: Nb of INTERNAL parameters of discretisation by U. */
3016 /* This is also the nb of root of Legendre polynom where discretization is done. */
3017 /* UROOTL: Table of parameters of discretisation ON (-1,1) by U.
3019 /* NBPNTV: Nb of INTERNAL parameters of discretisation by V. */
3020 /* This is also the nb of root of Legendre polynom where discretization is done. */
3021 /* IORDRV: Order of constraint imposed at the extremities of iso-V */
3022 /* = 0, calculate the extremities of iso-V */
3023 /* = 1, calculate, additionally, the 1st derivative in the direction of iso-V */
3024 /* = 2, calculate, additionally, the 2nd derivative in the direction of iso-V */
3025 /* SOTBU1: Table of NBPNTU/2 sums of 2 index points */
3026 /* NBPNTU-II+1 and II, for II = 1, NBPNTU/2 on iso-V0. */
3027 /* SOTBU2: Table of NBPNTV/2 sums of 2 index points */
3028 /* NBPNTU-II+1 and II, for II = 1, NBPNTU/2 on iso-V1. */
3029 /* DITBU1: Table of NBPNTU/2 differences of 2 index points */
3030 /* NBPNTU-II+1 and II, for II = 1, NBPNTU/2 on iso-V0. */
3031 /* DITBU2: Table of NBPNTU/2 differences of 2 index points */
3032 /* NBPNTU-II+1 and II, for II = 1, NBPNTU/2 on iso-V1. */
3033 /* SOSOTB: Preinitialized table (input/output argument). */
3034 /* DISOTB: Preinitialized table (input/output argument). */
3035 /* SODITB: Preinitialized table (input/output argument). */
3036 /* DIDITB: Preinitialized table (input/output argument) */
3038 /* OUTPUT ARGUMENTS : */
3039 /* ------------------- */
3040 /* FPNTAB: Auxiliary table. */
3041 /* UHERMT: Table of 2*(IORDRU+1) coeff. of 2*(IORDRU+1) polynoms of Hermite. */
3042 /* SOSOTB: Table where the terms of constraints are added */
3043 /* C(ui,vj) + C(ui,-vj) + C(-ui,vj) + C(-ui,-vj) */
3044 /* with ui and vj positive roots of the Legendre polynom */
3045 /* of degree NBPNTU and NBPNTV respectively. */
3046 /* DISOTB: Table where the terms of constraints are added */
3047 /* C(ui,vj) + C(ui,-vj) - C(-ui,vj) - C(-ui,-vj) */
3048 /* with ui and vj positive roots of the polynom of Legendre */
3049 /* of degree NBPNTU and NBPNTV respectively. */
3050 /* SODITB: Table where the terms of constraints are added */
3051 /* C(ui,vj) - C(ui,-vj) + C(-ui,vj) - C(-ui,-vj) */
3052 /* with ui and vj positive roots of the polynom of Legendre */
3053 /* of degree NBPNTU and NBPNTV respectively. */
3054 /* DIDITB: Table where the terms of constraints are added */
3055 /* C(ui,vj) - C(ui,-vj) - C(-ui,vj) + C(-ui,-vj) */
3056 /* with ui and vj positive roots of the polynom of Legendre */
3057 /* of degree NBPNTU and NBPNTV respectively. */
3059 /* COMMONS USED : */
3060 /* ---------------- */
3062 /* REFERENCES CALLED : */
3063 /* ----------------------- */
3065 /* DESCRIPTION/NOTES/LIMITATIONS : */
3066 /* ----------------------------------- */
3068 /* $ HISTORIQUE DES MODIFICATIONS : */
3069 /* -------------------------------- */
3070 /* 08-08-1991: RBD; Creation. */
3072 /* **********************************************************************
3075 /* Name of the routine */
3078 /* Parameter adjustments */
3080 diditb_dim1 = *nbpntu / 2 + 1;
3081 diditb_dim2 = *nbpntv / 2 + 1;
3082 diditb_offset = diditb_dim1 * diditb_dim2;
3083 diditb -= diditb_offset;
3084 disotb_dim1 = *nbpntu / 2;
3085 disotb_dim2 = *nbpntv / 2;
3086 disotb_offset = disotb_dim1 * (disotb_dim2 + 1) + 1;
3087 disotb -= disotb_offset;
3088 soditb_dim1 = *nbpntu / 2;
3089 soditb_dim2 = *nbpntv / 2;
3090 soditb_offset = soditb_dim1 * (soditb_dim2 + 1) + 1;
3091 soditb -= soditb_offset;
3092 sosotb_dim1 = *nbpntu / 2 + 1;
3093 sosotb_dim2 = *nbpntv / 2 + 1;
3094 sosotb_offset = sosotb_dim1 * sosotb_dim2;
3095 sosotb -= sosotb_offset;
3096 uhermt_dim1 = (*iordru << 1) + 2;
3097 uhermt_offset = uhermt_dim1;
3098 uhermt -= uhermt_offset;
3099 fpntab_dim1 = *nbpntu;
3100 fpntab_offset = fpntab_dim1 + 1;
3101 fpntab -= fpntab_offset;
3102 ditbu2_dim1 = *nbpntv / 2 + 1;
3103 ditbu2_dim2 = *ndimen;
3104 ditbu2_offset = ditbu2_dim1 * (ditbu2_dim2 + 1);
3105 ditbu2 -= ditbu2_offset;
3106 ditbu1_dim1 = *nbpntv / 2 + 1;
3107 ditbu1_dim2 = *ndimen;
3108 ditbu1_offset = ditbu1_dim1 * (ditbu1_dim2 + 1);
3109 ditbu1 -= ditbu1_offset;
3110 sotbu2_dim1 = *nbpntv / 2 + 1;
3111 sotbu2_dim2 = *ndimen;
3112 sotbu2_offset = sotbu2_dim1 * (sotbu2_dim2 + 1);
3113 sotbu2 -= sotbu2_offset;
3114 sotbu1_dim1 = *nbpntv / 2 + 1;
3115 sotbu1_dim2 = *ndimen;
3116 sotbu1_offset = sotbu1_dim1 * (sotbu1_dim2 + 1);
3117 sotbu1 -= sotbu1_offset;
3120 ibb = AdvApp2Var_SysBase::mnfndeb_();
3122 AdvApp2Var_SysBase::mgenmsg_("MMA2CD3", 7L);
3125 /* ------------------- Discretization of polynoms of Hermit -----------
3128 ncfhu = (*iordru + 1) << 1;
3130 for (ii = 1; ii <= i__1; ++ii) {
3132 for (kk = 1; kk <= i__2; ++kk) {
3133 AdvApp2Var_MathBase::mmmpocur_(&ncfhu,
3136 &uhermt[ii * uhermt_dim1],
3138 &fpntab[kk + ii * fpntab_dim1]);
3144 /* ---- The discretizations of polynoms of constraints are subtracted ----
3147 nvroo = *nbpntv / 2;
3148 nuroo = *nbpntu / 2;
3151 for (nd = 1; nd <= i__1; ++nd) {
3153 for (ii = 1; ii <= i__2; ++ii) {
3156 for (jj = 1; jj <= i__3; ++jj) {
3157 bid1 = sotbu1[jj + (nd + ii * sotbu1_dim2) * sotbu1_dim1];
3158 bid2 = sotbu2[jj + (nd + ii * sotbu2_dim2) * sotbu2_dim1];
3159 bid3 = ditbu1[jj + (nd + ii * ditbu1_dim2) * ditbu1_dim1];
3160 bid4 = ditbu2[jj + (nd + ii * ditbu2_dim2) * ditbu2_dim1];
3162 for (kk = 1; kk <= i__4; ++kk) {
3163 kkp = (*nbpntu + 1) / 2 + kk;
3164 kkm = nuroo - kk + 1;
3165 sosotb[kk + (jj + nd * sosotb_dim2) * sosotb_dim1] =
3166 sosotb[kk + (jj + nd * sosotb_dim2) * sosotb_dim1]
3167 - bid1 * (fpntab[kkp + ((ii << 1) - 1) *
3168 fpntab_dim1] + fpntab[kkm + ((ii << 1) - 1) *
3169 fpntab_dim1]) - bid2 * (fpntab[kkp + (ii << 1) *
3170 fpntab_dim1] + fpntab[kkm + (ii << 1) *
3172 disotb[kk + (jj + nd * disotb_dim2) * disotb_dim1] =
3173 disotb[kk + (jj + nd * disotb_dim2) * disotb_dim1]
3174 - bid1 * (fpntab[kkp + ((ii << 1) - 1) *
3175 fpntab_dim1] - fpntab[kkm + ((ii << 1) - 1) *
3176 fpntab_dim1]) - bid2 * (fpntab[kkp + (ii << 1) *
3177 fpntab_dim1] - fpntab[kkm + (ii << 1) *
3179 soditb[kk + (jj + nd * soditb_dim2) * soditb_dim1] =
3180 soditb[kk + (jj + nd * soditb_dim2) * soditb_dim1]
3181 - bid3 * (fpntab[kkp + ((ii << 1) - 1) *
3182 fpntab_dim1] + fpntab[kkm + ((ii << 1) - 1) *
3183 fpntab_dim1]) - bid4 * (fpntab[kkp + (ii << 1) *
3184 fpntab_dim1] + fpntab[kkm + (ii << 1) *
3186 diditb[kk + (jj + nd * diditb_dim2) * diditb_dim1] =
3187 diditb[kk + (jj + nd * diditb_dim2) * diditb_dim1]
3188 - bid3 * (fpntab[kkp + ((ii << 1) - 1) *
3189 fpntab_dim1] - fpntab[kkm + ((ii << 1) - 1) *
3190 fpntab_dim1]) - bid4 * (fpntab[kkp + (ii << 1) *
3191 fpntab_dim1] - fpntab[kkm + (ii << 1) *
3200 /* ------------ Case when the discretization is done only on the roots */
3201 /* ---------- of Legendre polynom of uneven degree, 0 is root */
3205 if (*nbpntu % 2 == 1) {
3207 for (ii = 1; ii <= i__2; ++ii) {
3209 for (jj = 1; jj <= i__3; ++jj) {
3210 bid1 = sotbu1[jj + (nd + ii * sotbu1_dim2) * sotbu1_dim1]
3211 * fpntab[nuroo + 1 + ((ii << 1) - 1) *
3212 fpntab_dim1] + sotbu2[jj + (nd + ii * sotbu2_dim2)
3213 * sotbu2_dim1] * fpntab[nuroo + 1 + (ii << 1) *
3215 sosotb[(jj + nd * sosotb_dim2) * sosotb_dim1] -= bid1;
3216 bid2 = ditbu1[jj + (nd + ii * ditbu1_dim2) * ditbu1_dim1]
3217 * fpntab[nuroo + 1 + ((ii << 1) - 1) *
3218 fpntab_dim1] + ditbu2[jj + (nd + ii * ditbu2_dim2)
3219 * ditbu2_dim1] * fpntab[nuroo + 1 + (ii << 1) *
3221 diditb[(jj + nd * diditb_dim2) * diditb_dim1] -= bid2;
3228 if (*nbpntv % 2 == 1) {
3230 for (ii = 1; ii <= i__2; ++ii) {
3232 for (kk = 1; kk <= i__3; ++kk) {
3233 kkp = (*nbpntu + 1) / 2 + kk;
3234 kkm = nuroo - kk + 1;
3235 bid1 = sotbu1[(nd + ii * sotbu1_dim2) * sotbu1_dim1] * (
3236 fpntab[kkp + ((ii << 1) - 1) * fpntab_dim1] +
3237 fpntab[kkm + ((ii << 1) - 1) * fpntab_dim1]) +
3238 sotbu2[(nd + ii * sotbu2_dim2) * sotbu2_dim1] * (
3239 fpntab[kkp + (ii << 1) * fpntab_dim1] + fpntab[
3240 kkm + (ii << 1) * fpntab_dim1]);
3241 sosotb[kk + nd * sosotb_dim2 * sosotb_dim1] -= bid1;
3242 bid2 = sotbu1[(nd + ii * sotbu1_dim2) * sotbu1_dim1] * (
3243 fpntab[kkp + ((ii << 1) - 1) * fpntab_dim1] -
3244 fpntab[kkm + ((ii << 1) - 1) * fpntab_dim1]) +
3245 sotbu2[(nd + ii * sotbu2_dim2) * sotbu2_dim1] * (
3246 fpntab[kkp + (ii << 1) * fpntab_dim1] - fpntab[
3247 kkm + (ii << 1) * fpntab_dim1]);
3248 diditb[kk + nd * diditb_dim2 * diditb_dim1] -= bid2;
3255 if (*nbpntu % 2 == 1 && *nbpntv % 2 == 1) {
3257 for (ii = 1; ii <= i__2; ++ii) {
3258 bid1 = sotbu1[(nd + ii * sotbu1_dim2) * sotbu1_dim1] * fpntab[
3259 nuroo + 1 + ((ii << 1) - 1) * fpntab_dim1] + sotbu2[(
3260 nd + ii * sotbu2_dim2) * sotbu2_dim1] * fpntab[nuroo
3261 + 1 + (ii << 1) * fpntab_dim1];
3262 sosotb[nd * sosotb_dim2 * sosotb_dim1] -= bid1;
3271 /* ------------------------------ The End -------------------------------
3276 AdvApp2Var_SysBase::mgsomsg_("MMA2CD3", 7L);
3281 //=======================================================================
3282 //function : mma2cdi_
3284 //=======================================================================
3285 int AdvApp2Var_ApproxF2var::mma2cdi_( integer *ndimen,
3313 /* System generated locals */
3314 integer contr1_dim1, contr1_dim2, contr1_offset, contr2_dim1, contr2_dim2,
3315 contr2_offset, contr3_dim1, contr3_dim2, contr3_offset,
3316 contr4_dim1, contr4_dim2, contr4_offset, sosotb_dim1, sosotb_dim2,
3317 sosotb_offset, diditb_dim1, diditb_dim2, diditb_offset,
3318 soditb_dim1, soditb_dim2, soditb_offset, disotb_dim1, disotb_dim2,
3321 /* Local variables */
3324 doublereal* wrkar = 0;
3326 integer ibb, ier = 0;
3327 integer isz1, isz2, isz3, isz4;
3328 intptr_t ipt1, ipt2, ipt3, ipt4;
3333 /* **********************************************************************
3338 /* Discretisation on the parameters of polynomes of interpolation */
3339 /* of constraints of order IORDRE. */
3343 /* TOUS, AB_SPECIFI::CONTRAINTE&, DISCRETISATION, &POINT */
3345 //* INPUT ARGUMENTS : */
3346 /* ------------------ */
3347 /* NDIMEN: Dimension of the space. */
3348 /* NBPNTU: Nb of INTERNAL parameters of discretisation by U. */
3349 /* This is also the nb of root of Legendre polynom where discretization is done. */
3350 /* UROOTL: Table of parameters of discretisation ON (-1,1) by U.
3352 /* NBPNTV: Nb of INTERNAL parameters of discretisation by V. */
3353 /* This is also the nb of root of Legendre polynom where discretization is done. */
3354 /* VROOTL: Table of parameters of discretisation ON (-1,1) by V.
3356 /* IORDRV: Order of constraint imposed at the extremities of iso-U */
3357 /* = 0, calculate the extremities of iso-U */
3358 /* = 1, calculate, additionally, the 1st derivative in the direction of iso-U */
3359 /* = 2, calculate, additionally, the 2nd derivative in the direction of iso-U */
3360 /* IORDRU: Order of constraint imposed at the extremities of iso-V */
3361 /* = 0, calculate the extremities of iso-V */
3362 /* = 1, calculate, additionally, the 1st derivative in the direction of iso-V */
3363 /* = 2, calculate, additionally, the 2nd derivative in the direction of iso-V */
3364 /* CONTR1: Contains, if IORDRU and IORDRV>=0, the values at the */
3365 /* extremities of F(U0,V0) and its derivatives. */
3366 /* CONTR2: Contains, if IORDRU and IORDRV>=0, the values at the */
3367 /* extremities of F(U1,V0) and its derivatives. */
3368 /* CONTR3: Contains, if IORDRU and IORDRV>=0, the values at the */
3369 /* extremities of F(U0,V1) and its derivatives. */
3370 /* CONTR4: Contains, if IORDRU and IORDRV>=0, the values at the */
3371 /* extremities of F(U1,V1) and its derivatives. */
3372 /* SOTBU1: Table of NBPNTU/2 sums of 2 index points */
3373 /* NBPNTU-II+1 and II, for II = 1, NBPNTU/2 on iso-V0. */
3374 /* SOTBU2: Table of NBPNTV/2 sums of 2 index points */
3375 /* NBPNTU-II+1 and II, for II = 1, NBPNTU/2 on iso-V1. */
3376 /* DITBU1: Table of NBPNTU/2 differences of 2 index points */
3377 /* NBPNTU-II+1 and II, for II = 1, NBPNTU/2 on iso-V0. */
3378 /* DITBU2: Table of NBPNTU/2 differences of 2 index points */
3379 /* NBPNTU-II+1 and II, for II = 1, NBPNTU/2 on iso-V1. */
3380 /* SOTBV1: Table of NBPNTV/2 sums of 2 index points */
3381 /* NBPNTV-II+1 and II, for II = 1, NBPNTV/2 on iso-V0. */
3382 /* SOTBV2: Table of NBPNTV/2 sums of 2 index points */
3383 /* NBPNTV-II+1 and II, for II = 1, NBPNTV/2 on iso-V1. */
3384 /* DITBV1: Table of NBPNTV/2 differences of 2 index points */
3385 /* NBPNTV-II+1 and II, for II = 1, NBPNTV/2 on iso-V0. */
3386 /* DITBV2: Table of NBPNTV/2 differences of 2 index points */
3387 /* NBPNTV-II+1 and II, for II = 1, NBPNTV/2 on iso-V1. */
3388 /* SOSOTB: Preinitialized table (input/output argument). */
3389 /* DISOTB: Preinitialized table (input/output argument). */
3390 /* SODITB: Preinitialized table (input/output argument). */
3391 /* DIDITB: Preinitialized table (input/output argument) */
3393 /* ARGUMENTS DE SORTIE : */
3394 /* ------------------- */
3395 /* SOSOTB: Table where the terms of constraints are added */
3396 /* C(ui,vj) + C(ui,-vj) + C(-ui,vj) + C(-ui,-vj) */
3397 /* with ui and vj positive roots of the Legendre polynom */
3398 /* of degree NBPNTU and NBPNTV respectively. */
3399 /* DISOTB: Table where the terms of constraints are added */
3400 /* C(ui,vj) + C(ui,-vj) - C(-ui,vj) - C(-ui,-vj) */
3401 /* with ui and vj positive roots of the polynom of Legendre */
3402 /* of degree NBPNTU and NBPNTV respectively. */
3403 /* SODITB: Table where the terms of constraints are added */
3404 /* C(ui,vj) - C(ui,-vj) + C(-ui,vj) - C(-ui,-vj) */
3405 /* with ui and vj positive roots of the polynom of Legendre */
3406 /* of degree NBPNTU and NBPNTV respectively. */
3407 /* DIDITB: Table where the terms of constraints are added */
3408 /* C(ui,vj) - C(ui,-vj) - C(-ui,vj) + C(-ui,-vj) */
3409 /* with ui and vj positive roots of the polynom of Legendre */
3410 /* of degree NBPNTU and NBPNTV respectively. */
3411 /* IERCOD: = 0, OK, */
3412 /* = 1, Value or IORDRV or IORDRU is out of allowed values. */
3413 /* =13, Pb of dynamic allocation. */
3415 /* COMMONS USED : */
3416 /* ---------------- */
3418 /* REFERENCES CALLED : */
3419 /* -------------------- */
3421 /* DESCRIPTION/NOTES/LIMITATIONS : */
3422 /* ------------------------------- */
3425 /* **********************************************************************
3428 /* The name of the routine */
3431 /* Parameter adjustments */
3433 diditb_dim1 = *nbpntu / 2 + 1;
3434 diditb_dim2 = *nbpntv / 2 + 1;
3435 diditb_offset = diditb_dim1 * diditb_dim2;
3436 diditb -= diditb_offset;
3437 disotb_dim1 = *nbpntu / 2;
3438 disotb_dim2 = *nbpntv / 2;
3439 disotb_offset = disotb_dim1 * (disotb_dim2 + 1) + 1;
3440 disotb -= disotb_offset;
3441 soditb_dim1 = *nbpntu / 2;
3442 soditb_dim2 = *nbpntv / 2;
3443 soditb_offset = soditb_dim1 * (soditb_dim2 + 1) + 1;
3444 soditb -= soditb_offset;
3445 sosotb_dim1 = *nbpntu / 2 + 1;
3446 sosotb_dim2 = *nbpntv / 2 + 1;
3447 sosotb_offset = sosotb_dim1 * sosotb_dim2;
3448 sosotb -= sosotb_offset;
3450 contr4_dim1 = *ndimen;
3451 contr4_dim2 = *iordru + 2;
3452 contr4_offset = contr4_dim1 * (contr4_dim2 + 1) + 1;
3453 contr4 -= contr4_offset;
3454 contr3_dim1 = *ndimen;
3455 contr3_dim2 = *iordru + 2;
3456 contr3_offset = contr3_dim1 * (contr3_dim2 + 1) + 1;
3457 contr3 -= contr3_offset;
3458 contr2_dim1 = *ndimen;
3459 contr2_dim2 = *iordru + 2;
3460 contr2_offset = contr2_dim1 * (contr2_dim2 + 1) + 1;
3461 contr2 -= contr2_offset;
3462 contr1_dim1 = *ndimen;
3463 contr1_dim2 = *iordru + 2;
3464 contr1_offset = contr1_dim1 * (contr1_dim2 + 1) + 1;
3465 contr1 -= contr1_offset;
3474 AdvApp2Var_SysBase anAdvApp2Var_SysBase;
3477 ibb = AdvApp2Var_SysBase::mnfndeb_();
3479 AdvApp2Var_SysBase::mgenmsg_("MMA2CDI", 7L);
3483 if (*iordru < -1 || *iordru > 2) {
3486 if (*iordrv < -1 || *iordrv > 2) {
3490 /* ------------------------- Set to zero --------------------------------
3493 ilong = (*nbpntu / 2 + 1) * (*nbpntv / 2 + 1) * *ndimen;
3494 AdvApp2Var_SysBase::mvriraz_(&ilong, &sosotb[sosotb_offset]);
3495 AdvApp2Var_SysBase::mvriraz_(&ilong, &diditb[diditb_offset]);
3496 ilong = *nbpntu / 2 * (*nbpntv / 2) * *ndimen;
3497 AdvApp2Var_SysBase::mvriraz_(&ilong, &soditb[soditb_offset]);
3498 AdvApp2Var_SysBase::mvriraz_(&ilong, &disotb[disotb_offset]);
3499 if (*iordru == -1 && *iordrv == -1) {
3505 isz1 = ((*iordru + 1) << 2) * (*iordru + 1);
3506 isz2 = ((*iordrv + 1) << 2) * (*iordrv + 1);
3507 isz3 = ((*iordru + 1) << 1) * *nbpntu;
3508 isz4 = ((*iordrv + 1) << 1) * *nbpntv;
3509 iszwr = isz1 + isz2 + isz3 + isz4;
3510 anAdvApp2Var_SysBase.mcrrqst_(&c__8, &iszwr, wrkar, &iofwr, &ier);
3519 if (*iordru >= 0 && *iordru <= 2) {
3521 /* --- Return 2*(IORDRU+1) coeff of 2*(IORDRU+1) polynoms of Hermite
3524 AdvApp2Var_ApproxF2var::mma1her_(iordru, &wrkar[ipt1], iercod);
3529 /* ---- Subract discretizations of polynoms of constraints
3532 mma2cd3_(ndimen, nbpntu, &urootl[1], nbpntv, iordru, &sotbu1[1], &
3533 sotbu2[1], &ditbu1[1], &ditbu2[1], &wrkar[ipt3], &wrkar[ipt1],
3534 &sosotb[sosotb_offset], &soditb[soditb_offset], &disotb[
3535 disotb_offset], &diditb[diditb_offset]);
3538 if (*iordrv >= 0 && *iordrv <= 2) {
3540 /* --- Return 2*(IORDRV+1) coeff of 2*(IORDRV+1) polynoms of Hermite
3543 AdvApp2Var_ApproxF2var::mma1her_(iordrv, &wrkar[ipt2], iercod);
3548 /* ---- Subtract discretisations of polynoms of constraint
3551 mma2cd2_(ndimen, nbpntu, nbpntv, &vrootl[1], iordrv, &sotbv1[1], &
3552 sotbv2[1], &ditbv1[1], &ditbv2[1], &wrkar[ipt4], &wrkar[ipt2],
3553 &sosotb[sosotb_offset], &soditb[soditb_offset], &disotb[
3554 disotb_offset], &diditb[diditb_offset]);
3557 /* --------------- Subtract constraints of corners ----------------
3560 if (*iordru >= 0 && *iordrv >= 0) {
3561 mma2cd1_(ndimen, nbpntu, &urootl[1], nbpntv, &vrootl[1], iordru,
3562 iordrv, &contr1[contr1_offset], &contr2[contr2_offset], &
3563 contr3[contr3_offset], &contr4[contr4_offset], &wrkar[ipt3], &
3564 wrkar[ipt4], &wrkar[ipt1], &wrkar[ipt2], &sosotb[
3565 sosotb_offset], &soditb[soditb_offset], &disotb[disotb_offset]
3566 , &diditb[diditb_offset]);
3570 /* ------------------------------ The End -------------------------------
3572 /* --> IORDRE is not within the autorised diapason. */
3576 /* --> PB of dynamic allocation. */
3583 anAdvApp2Var_SysBase.mcrdelt_(&c__8, &iszwr, wrkar, &iofwr, &ier);
3588 AdvApp2Var_SysBase::maermsg_("MMA2CDI", iercod, 7L);
3590 AdvApp2Var_SysBase::mgsomsg_("MMA2CDI", 7L);
3595 //=======================================================================
3596 //function : mma2ce1_
3598 //=======================================================================
3599 int AdvApp2Var_ApproxF2var::mma2ce1_(integer *numdec,
3629 /* System generated locals */
3630 integer sosotb_dim1, sosotb_dim2, sosotb_offset, disotb_dim1, disotb_dim2,
3631 disotb_offset, soditb_dim1, soditb_dim2, soditb_offset,
3632 diditb_dim1, diditb_dim2, diditb_offset, patjac_dim1, patjac_dim2,
3635 /* Local variables */
3638 doublereal* wrkar = 0;
3641 integer isz1, isz2, isz3, isz4, isz5, isz6, isz7;
3642 intptr_t ipt1, ipt2, ipt3, ipt4, ipt5, ipt6, ipt7;
3646 /* **********************************************************************
3651 /* Calculation of coefficients of polynomial approximation of degree */
3652 /* (NDJACU,NDJACV) of a function F(u,v), starting from its */
3653 /* discretization on roots of Legendre polynom of degree */
3654 /* NBPNTU by U and NBPNTV by V. */
3658 /* TOUS,AB_SPECIFI::FONCTION&,APPROXIMATION,&POLYNOME,&ERREUR */
3660 /* INPUT ARGUMENTS : */
3661 /* ------------------ */
3662 /* NUMDEC: Indicates if it is POSSIBLE to cut function F(u,v). */
3663 /* = 5, It is POSSIBLE to cut by U or by V or in both directions simultaneously. */
3664 /* = 4, It is POSSIBLE to cut by U or by V BUT NOT in both */
3665 /* directions simultaneously (cutting by V is preferable). */
3666 /* = 3, It is POSSIBLE to cut by U or by V BUT NOT in both */
3667 /* directions simultaneously (cutting by U is preferable). */
3668 /* = 2, It is POSSIBLE to cut only by V (i.e. insert parameter */
3669 /* of cutting Vj). */
3670 /* = 1, It is POSSIBLE to cut only by U (i.e. insert parameter */
3671 /* of cutting Ui). */
3672 /* = 0, It is not POSSIBLE to cut anything */
3673 /* NDIMEN: Dimension of the space. */
3674 /* NBSESP: Nb of independent sub-spaces on which the errors are calculated. */
3675 /* NDIMSE: Table of dimensions of each of sub-spaces. */
3676 /* NDMINU: Minimum degree by U to be preserved for the approximation. */
3677 /* NDMINV: Minimum degree by V to be preserved for the approximation. */
3678 /* NDGULI: Limit of nb of coefficients by U of the solution. */
3679 /* NDGVLI: Limit of nb of coefficients by V of the solution. */
3680 /* NDJACU: Max degree of the polynom of approximation by U. */
3681 /* The representation in the orthogonal base starts from degree */
3682 /* 0 to degree NDJACU-2*(IORDRU+1). The polynomial base is the base of */
3683 /* Jacobi of order -1 (Legendre), 0, 1 or 2. */
3684 /* It is required that 2*IORDRU+1 <= NDMINU <= NDGULI < NDJACU */
3685 /* NDJACV: Max degree of the polynom of approximation by V. */
3686 /* The representation in the orthogonal base starts from degree */
3687 /* 0 to degree NDJACV-2*(IORDRV+1). The polynomial base is */
3688 /* the base of Jacobi of order -1 (Legendre), 0, 1 or 2 */
3689 /* It is required that 2*IORDRV+1 <= NDMINV <= NDGVLI < NDJACV */
3690 /* IORDRU: Order of the Jacobi base (-1,0,1 or 2) by U. Corresponds */
3691 /* to the step of constraints C0, C1 or C2. */
3692 /* IORDRV: Order of the Jacobi base (-1,0,1 or 2) by U. Corresponds */
3693 /* to the step of constraints C0, C1 or C2. */
3694 /* NBPNTU: Degree of Legendre polynom on the roots which of are */
3695 /* calculated the coefficients of integration by u */
3696 /* by Gauss method. It is required that NBPNTU = 30, 40, */
3697 /* 50 or 61 and NDJACU-2*(IORDRU+1) < NBPNTU. */
3698 /* NBPNTV: Degree of Legendre polynom on the roots which of are */
3699 /* calculated the coefficients of integration by u */
3700 /* by Gauss method. It is required that NBPNTV = 30, 40, */
3701 /* 50 or 61 and NDJACV-2*(IORDRV+1) < NBPNTV. */
3702 /* EPSAPR: Table of NBSESP tolerances imposed on each sub-spaces. */
3703 /* SOSOTB: Table of F(ui,vj) + F(ui,-vj) + F(-ui,vj) + F(-ui,-vj) */
3704 /* with ui and vj - positive roots of the Legendre polynom */
3705 /* of degree NBPNTU and NBPNTV respectively. Additionally, */
3706 /* table SOSOTB(0,j) contains F(0,vj) + F(0,-vj), */
3707 /* table SOSOTB(i,0) contains F(ui,0) + F(-ui,0) and */
3708 /* SOSOTB(0,0) contains F(0,0). */
3709 /* DISOTB: Table of F(ui,vj) + F(ui,-vj) - F(-ui,vj) - F(-ui,-vj) */
3710 /* with ui and vj positive roots of Legendre polynom */
3711 /* of degree NBPNTU and NBPNTV respectively. */
3712 /* SODITB: Table of F(ui,vj) - F(ui,-vj) + F(-ui,vj) - F(-ui,-vj) */
3713 /* with ui and vj positive roots of Legendre polynom */
3714 /* of degree NBPNTU and NBPNTV respectively. */
3715 /* DIDITB: Table of F(ui,vj) - F(ui,-vj) - F(-ui,vj) + F(-ui,-vj) */
3716 /* with ui and vj positive roots of Legendre polynom */
3717 /* of degree NBPNTU and NBPNTV respectively. Additionally, */
3718 /* table DIDITB(0,j) contains F(0,vj) - F(0,-vj), */
3719 /* and table DIDITB(i,0) contains F(ui,0) - F(-ui,0). */
3721 /* OUTPUT ARGUMENTS */
3722 /* --------------- */
3723 /* PATJAC: Table of coefficients of polynom P(u,v) of approximation */
3724 /* of F(u,v) with eventually taking into account of */
3725 /* constraints. P(u,v) is of degree (NDJACU,NDJACV). */
3726 /* This table contains other coeff if ITYDEC = 0. */
3727 /* ERRMAX: For 1<=i<=NBSESP, ERRMAX(i) contains max errors */
3728 /* on each of sub-spaces SI ITYDEC = 0. */
3729 /* ERRMOY: Contains average errors for each of NBSESP sub-spaces SI ITYDEC = 0. */
3730 /* NDEGPU: Degree by U for square PATJAC. Valable if ITYDEC=0. */
3731 /* NDEGPV: Degree by V for square PATJAC. Valable if ITYDEC=0. */
3732 /* ITYDEC: Shows if it is NECESSARY to cut again function F(u,v). */
3733 /* = 0, it is not NECESSARY to cut anything, PATJAC is OK. */
3734 /* = 1, it is NECESSARY to cut only by U (i.e. insert parameter of cutting Ui). */
3735 /* = 2, it is NECESSARY to cut only by V (i.e. insert parameter of cutting Vj). */
3736 /* = 3, it is NECESSARY to cut both by U AND by V. */
3737 /* IERCOD: Error code. */
3738 /* = 0, Everything is OK. */
3739 /* = -1, There is the best possible solution, but the */
3740 /* user tolerance is not satisfactory (3*only) */
3741 /* = 1, Incoherent entries. */
3743 /* COMMONS USED : */
3744 /* ---------------- */
3746 /* REFERENCES CALLED : */
3747 /* --------------------- */
3749 /* DESCRIPTION/NOTES/LIMITATIONS : */
3750 /* ------------------------------- */
3753 /* **********************************************************************
3755 /* Name of the routine */
3758 /* --------------------------- Initialisations --------------------------
3761 /* Parameter adjustments */
3766 patjac_dim1 = *ndjacu + 1;
3767 patjac_dim2 = *ndjacv + 1;
3768 patjac_offset = patjac_dim1 * patjac_dim2;
3769 patjac -= patjac_offset;
3770 diditb_dim1 = *nbpntu / 2 + 1;
3771 diditb_dim2 = *nbpntv / 2 + 1;
3772 diditb_offset = diditb_dim1 * diditb_dim2;
3773 diditb -= diditb_offset;
3774 soditb_dim1 = *nbpntu / 2;
3775 soditb_dim2 = *nbpntv / 2;
3776 soditb_offset = soditb_dim1 * (soditb_dim2 + 1) + 1;
3777 soditb -= soditb_offset;
3778 disotb_dim1 = *nbpntu / 2;
3779 disotb_dim2 = *nbpntv / 2;
3780 disotb_offset = disotb_dim1 * (disotb_dim2 + 1) + 1;
3781 disotb -= disotb_offset;
3782 sosotb_dim1 = *nbpntu / 2 + 1;
3783 sosotb_dim2 = *nbpntv / 2 + 1;
3784 sosotb_offset = sosotb_dim1 * sosotb_dim2;
3785 sosotb -= sosotb_offset;
3788 ldbg = AdvApp2Var_SysBase::mnfndeb_() >= 3;
3790 AdvApp2Var_SysBase::mgenmsg_("MMA2CE1", 7L);
3795 isz1 = (*nbpntu / 2 + 1) * (*ndjacu - ((*iordru + 1) << 1) + 1);
3796 isz2 = (*nbpntv / 2 + 1) * (*ndjacv - ((*iordrv + 1) << 1) + 1);
3797 isz3 = (*nbpntv / 2 + 1) * (*ndjacu - ((*iordru + 1) << 1) + 1) * *ndimen;
3798 isz4 = *nbpntv / 2 * (*ndjacu - ((*iordru + 1) << 1) + 1) * *ndimen;
3799 isz5 = *ndjacu + 1 - ((*iordru + 1) << 1);
3800 isz6 = *ndjacv + 1 - ((*iordrv + 1) << 1);
3801 isz7 = *ndimen << 2;
3802 iszwr = isz1 + isz2 + isz3 + isz4 + isz5 + isz6 + isz7;
3803 AdvApp2Var_SysBase anAdvApp2Var_SysBase;
3804 anAdvApp2Var_SysBase.mcrrqst_(&c__8, &iszwr, wrkar, &iofwr, &ier);
3816 /* ----------------- Return Gauss coefficients of integration ----------------
3819 AdvApp2Var_ApproxF2var::mmapptt_(ndjacu, nbpntu, iordru, &wrkar[ipt1], iercod);
3823 AdvApp2Var_ApproxF2var::mmapptt_(ndjacv, nbpntv, iordrv, &wrkar[ipt2], iercod);
3828 /* ------------------- Return max polynoms of Jacobi ------------
3831 AdvApp2Var_ApproxF2var::mma2jmx_(ndjacu, iordru, &wrkar[ipt5]);
3832 AdvApp2Var_ApproxF2var::mma2jmx_(ndjacv, iordrv, &wrkar[ipt6]);
3834 /* ------ Calculate the coefficients and their contribution to the error ----
3837 mma2ce2_(numdec, ndimen, nbsesp, &ndimse[1], ndminu, ndminv, ndguli,
3838 ndgvli, ndjacu, ndjacv, iordru, iordrv, nbpntu, nbpntv, &epsapr[1]
3839 , &sosotb[sosotb_offset], &disotb[disotb_offset], &soditb[
3840 soditb_offset], &diditb[diditb_offset], &wrkar[ipt1], &wrkar[ipt2]
3841 , &wrkar[ipt5], &wrkar[ipt6], &wrkar[ipt7], &wrkar[ipt3], &wrkar[
3842 ipt4], &patjac[patjac_offset], &errmax[1], &errmoy[1], ndegpu,
3843 ndegpv, itydec, iercod);
3849 /* ------------------------------ The end -------------------------------
3858 anAdvApp2Var_SysBase.mcrdelt_(&c__8, &iszwr, wrkar, &iofwr, &ier);
3863 AdvApp2Var_SysBase::maermsg_("MMA2CE1", iercod, 7L);
3865 AdvApp2Var_SysBase::mgsomsg_("MMA2CE1", 7L);
3870 //=======================================================================
3871 //function : mma2ce2_
3873 //=======================================================================
3874 int mma2ce2_(integer *numdec,
3909 /* System generated locals */
3910 integer sosotb_dim1, sosotb_dim2, sosotb_offset, disotb_dim1, disotb_dim2,
3911 disotb_offset, soditb_dim1, soditb_dim2, soditb_offset,
3912 diditb_dim1, diditb_dim2, diditb_offset, gssutb_dim1, gssvtb_dim1,
3913 chpair_dim1, chpair_dim2, chpair_offset, chimpr_dim1,
3914 chimpr_dim2, chimpr_offset, patjac_dim1, patjac_dim2,
3915 patjac_offset, vecerr_dim1, vecerr_offset, i__1, i__2, i__3, i__4;
3917 /* Local variables */
3919 integer idim, igsu, minu, minv, maxu, maxv, igsv;
3921 integer i2rdu, i2rdv, ndses, nd, ii, jj, kk, nu, nv;
3925 /* **********************************************************************
3929 /* Calculation of coefficients of polynomial approximation of degree */
3930 /* (NDJACU,NDJACV) of a function F(u,v), starting from its */
3931 /* discretization on roots of Legendre polynom of degree */
3932 /* NBPNTU by U and NBPNTV by V. */
3936 /* TOUS,AB_SPECIFI::FONCTION&,APPROXIMATION,&COEFFICIENT,&POLYNOME */
3938 /* INPUT ARGUMENTS : */
3939 /* ------------------ */
3940 /* NUMDEC: Indicates if it is POSSIBLE to cut function F(u,v). */
3941 /* = 5, It is POSSIBLE to cut by U or by V or in both directions simultaneously. */
3942 /* = 4, It is POSSIBLE to cut by U or by V BUT NOT in both */
3943 /* directions simultaneously (cutting by V is preferable). */
3944 /* = 3, It is POSSIBLE to cut by U or by V BUT NOT in both */
3945 /* directions simultaneously (cutting by U is preferable). */
3946 /* = 2, It is POSSIBLE to cut only by V (i.e. insert parameter */
3947 /* of cutting Vj). */
3948 /* = 1, It is POSSIBLE to cut only by U (i.e. insert parameter */
3949 /* of cutting Ui). */
3950 /* = 0, It is not POSSIBLE to cut anything */
3951 /* NDIMEN: Total dimension of the space. */
3952 /* NBSESP: Nb of independent sub-spaces on which the errors are calculated. */
3953 /* NDIMSE: Table of dimensions of each of sub-spaces. */
3954 /* NDMINU: Minimum degree by U to be preserved for the approximation. */
3955 /* NDMINV: Minimum degree by V to be preserved for the approximation. */
3956 /* NDGULI: Limit of nb of coefficients by U of the solution. */
3957 /* NDGVLI: Limit of nb of coefficients by V of the solution. */
3958 /* NDJACU: Max degree of the polynom of approximation by U. */
3959 /* The representation in the orthogonal base starts from degree */
3960 /* 0 to degree NDJACU-2*(IORDRU+1). The polynomial base is the base of */
3961 /* Jacobi of order -1 (Legendre), 0, 1 or 2. */
3962 /* It is required that 2*IORDRU+1 <= NDMINU <= NDGULI < NDJACU */
3963 /* NDJACV: Max degree of the polynom of approximation by V. */
3964 /* The representation in the orthogonal base starts from degree */
3965 /* 0 to degree NDJACV-2*(IORDRV+1). The polynomial base is */
3966 /* the base of Jacobi of order -1 (Legendre), 0, 1 or 2 */
3967 /* It is required that 2*IORDRV+1 <= NDMINV <= NDGVLI < NDJACV */
3968 /* IORDRU: Order of the Jacobi base (-1,0,1 or 2) by U. Corresponds */
3969 /* to the step of constraints C0, C1 or C2. */
3970 /* IORDRV: Order of the Jacobi base (-1,0,1 or 2) by U. Corresponds */
3971 /* to the step of constraints C0, C1 or C2. */
3972 /* NBPNTU: Degree of Legendre polynom on the roots which of are */
3973 /* calculated the coefficients of integration by u */
3974 /* by Gauss method. It is required that NBPNTU = 30, 40, */
3975 /* 50 or 61 and NDJACU-2*(IORDRU+1) < NBPNTU. */
3976 /* NBPNTV: Degree of Legendre polynom on the roots which of are */
3977 /* calculated the coefficients of integration by u */
3978 /* by Gauss method. It is required that NBPNTV = 30, 40, */
3979 /* 50 or 61 and NDJACV-2*(IORDRV+1) < NBPNTV. */
3980 /* EPSAPR: Table of NBSESP tolerances imposed on each sub-spaces. */
3981 /* SOSOTB: Table of F(ui,vj) + F(ui,-vj) + F(-ui,vj) + F(-ui,-vj) */
3982 /* with ui and vj - positive roots of the Legendre polynom */
3983 /* of degree NBPNTU and NBPNTV respectively. Additionally, */
3984 /* table SOSOTB(0,j) contains F(0,vj) + F(0,-vj), */
3985 /* table SOSOTB(i,0) contains F(ui,0) + F(-ui,0) and */
3986 /* SOSOTB(0,0) contains F(0,0). */
3987 /* DISOTB: Table of F(ui,vj) + F(ui,-vj) - F(-ui,vj) - F(-ui,-vj) */
3988 /* with ui and vj positive roots of Legendre polynom */
3989 /* of degree NBPNTU and NBPNTV respectively. */
3990 /* SODITB: Table of F(ui,vj) - F(ui,-vj) + F(-ui,vj) - F(-ui,-vj) */
3991 /* with ui and vj positive roots of Legendre polynom */
3992 /* of degree NBPNTU and NBPNTV respectively. */
3993 /* DIDITB: Table of F(ui,vj) - F(ui,-vj) - F(-ui,vj) + F(-ui,-vj) */
3994 /* with ui and vj positive roots of Legendre polynom */
3995 /* of degree NBPNTU and NBPNTV respectively. Additionally, */
3996 /* table DIDITB(0,j) contains F(0,vj) - F(0,-vj), */
3997 /* and table DIDITB(i,0) contains F(ui,0) - F(-ui,0). */
3998 /* GSSUTB: Table of coefficients of integration by Gauss method */
3999 /* by U: i varies from 0 to NBPNTU/2 and k varies from 0 to */
4000 /* NDJACU-2*(IORDRU+1). */
4001 /* GSSVTB: Table of coefficients of integration by Gauss method */
4002 /* by V: i varies from 0 to NBPNTV/2 and k varies from 0 to */
4003 /* NDJACV-2*(IORDRV+1). */
4004 /* XMAXJU: Maximum value of Jacobi polynoms of order IORDRU, */
4005 /* from degree 0 to degree NDJACU - 2*(IORDRU+1) */
4006 /* XMAXJV: Maximum value of Jacobi polynoms of order IORDRV, */
4007 /* from degree 0 to degree NDJACV - 2*(IORDRV+1) */
4009 /* OUTPUT ARGUMENTS : */
4010 /* ------------------- */
4011 /* VECERR: Auxiliary table. */
4012 /* CHPAIR: Auxiliary table of terms connected to degree NDJACU by U */
4013 /* to calculate the coeff. of approximation of EVEN degree by V. */
4014 /* CHIMPR: Auxiliary table of terms connected to degree NDJACU by U */
4015 /* to calculate the coeff. of approximation of UNEVEN degree by V. */
4016 /* PATJAC: Table of coefficients of polynom P(u,v) of approximation */
4017 /* of F(u,v) with eventually taking into account of */
4018 /* constraints. P(u,v) is of degree (NDJACU,NDJACV). */
4019 /* This table contains other coeff if ITYDEC = 0. */
4020 /* ERRMAX: For 1<=i<=NBSESP, ERRMAX(i) contains max errors */
4021 /* on each of sub-spaces SI ITYDEC = 0. */
4022 /* ERRMOY: Contains average errors for each of NBSESP sub-spaces SI ITYDEC = 0. */
4023 /* NDEGPU: Degree by U for square PATJAC. Valable if ITYDEC=0. */
4024 /* NDEGPV: Degree by V for square PATJAC. Valable if ITYDEC=0. */
4025 /* ITYDEC: Shows if it is NECESSARY to cut again function F(u,v). */
4026 /* = 0, it is not NECESSARY to cut anything, PATJAC is OK. */
4027 /* = 1, it is NECESSARY to cut only by U (i.e. insert parameter of cutting Ui). */
4028 /* = 2, it is NECESSARY to cut only by V (i.e. insert parameter of cutting Vj). */
4029 /* = 3, it is NECESSARY to cut both by U AND by V. */
4030 /* IERCOD: Error code. */
4031 /* = 0, Everything is OK. */
4032 /* = -1, There is the best possible solution, but the */
4033 /* user tolerance is not satisfactory (3*only) */
4034 /* = 1, Incoherent entries. */
4036 /* COMMONS USED : */
4037 /* ---------------- */
4039 /* REFERENCES CALLED : */
4040 /* --------------------- */
4042 /* DESCRIPTION/NOTES/LIMITATIONS : */
4044 /* **********************************************************************
4046 /* Name of the routine */
4049 /* --------------------------- Initialisations --------------------------
4052 /* Parameter adjustments */
4053 vecerr_dim1 = *ndimen;
4054 vecerr_offset = vecerr_dim1 + 1;
4055 vecerr -= vecerr_offset;
4060 patjac_dim1 = *ndjacu + 1;
4061 patjac_dim2 = *ndjacv + 1;
4062 patjac_offset = patjac_dim1 * patjac_dim2;
4063 patjac -= patjac_offset;
4064 gssutb_dim1 = *nbpntu / 2 + 1;
4065 chimpr_dim1 = *nbpntv / 2;
4066 chimpr_dim2 = *ndjacu - ((*iordru + 1) << 1) + 1;
4067 chimpr_offset = chimpr_dim1 * chimpr_dim2 + 1;
4068 chimpr -= chimpr_offset;
4069 chpair_dim1 = *nbpntv / 2 + 1;
4070 chpair_dim2 = *ndjacu - ((*iordru + 1) << 1) + 1;
4071 chpair_offset = chpair_dim1 * chpair_dim2;
4072 chpair -= chpair_offset;
4073 gssvtb_dim1 = *nbpntv / 2 + 1;
4074 diditb_dim1 = *nbpntu / 2 + 1;
4075 diditb_dim2 = *nbpntv / 2 + 1;
4076 diditb_offset = diditb_dim1 * diditb_dim2;
4077 diditb -= diditb_offset;
4078 soditb_dim1 = *nbpntu / 2;
4079 soditb_dim2 = *nbpntv / 2;
4080 soditb_offset = soditb_dim1 * (soditb_dim2 + 1) + 1;
4081 soditb -= soditb_offset;
4082 disotb_dim1 = *nbpntu / 2;
4083 disotb_dim2 = *nbpntv / 2;
4084 disotb_offset = disotb_dim1 * (disotb_dim2 + 1) + 1;
4085 disotb -= disotb_offset;
4086 sosotb_dim1 = *nbpntu / 2 + 1;
4087 sosotb_dim2 = *nbpntv / 2 + 1;
4088 sosotb_offset = sosotb_dim1 * sosotb_dim2;
4089 sosotb -= sosotb_offset;
4092 ldbg = AdvApp2Var_SysBase::mnfndeb_() >= 3;
4094 AdvApp2Var_SysBase::mgenmsg_("MMA2CE2", 7L);
4096 /* --> A priori everything is OK */
4098 /* --> test of inputs */
4099 if (*numdec < 0 || *numdec > 5) {
4102 if ((*iordru << 1) + 1 > *ndminu) {
4105 if (*ndminu > *ndguli) {
4108 if (*ndguli >= *ndjacu) {
4111 if ((*iordrv << 1) + 1 > *ndminv) {
4114 if (*ndminv > *ndgvli) {
4117 if (*ndgvli >= *ndjacv) {
4120 /* --> A priori, no cuts to be done */
4122 /* --> Min. degrees to return: NDMINU,NDMINV */
4125 /* --> For the moment, max errors are null */
4126 AdvApp2Var_SysBase::mvriraz_(nbsesp, &errmax[1]);
4128 AdvApp2Var_SysBase::mvriraz_(&nd, &vecerr[vecerr_offset]);
4129 /* --> and the square, too. */
4130 nd = (*ndjacu + 1) * (*ndjacv + 1) * *ndimen;
4131 AdvApp2Var_SysBase::mvriraz_(&nd, &patjac[patjac_offset]);
4133 i2rdu = (*iordru + 1) << 1;
4134 i2rdv = (*iordrv + 1) << 1;
4136 /* **********************************************************************
4138 /* -------------------- HERE IT IS POSSIBLE TO CUT ----------------------
4140 /* **********************************************************************
4143 if (*numdec > 0 && *numdec <= 5) {
4145 /* ******************************************************************
4147 /* ---------------------- Calculate coeff of zone 4 -------------
4161 /* ---------------- Calculate the terms connected to degree by U ---------
4165 for (nd = 1; nd <= i__1; ++nd) {
4167 for (kk = minu; kk <= i__2; ++kk) {
4169 mma2cfu_(&kk, nbpntu, nbpntv, &sosotb[nd * sosotb_dim2 *
4170 sosotb_dim1], &disotb[(nd * disotb_dim2 + 1) *
4171 disotb_dim1 + 1], &soditb[(nd * soditb_dim2 + 1) *
4172 soditb_dim1 + 1], &diditb[nd * diditb_dim2 *
4173 diditb_dim1], &gssutb[igsu * gssutb_dim1], &chpair[(
4174 igsu + nd * chpair_dim2) * chpair_dim1], &chimpr[(
4175 igsu + nd * chimpr_dim2) * chimpr_dim1 + 1]);
4181 /* ------------------- Calculate the coefficients of PATJAC ------------
4184 igsu = minu - i2rdu;
4186 for (jj = minv; jj <= i__1; ++jj) {
4189 for (nd = 1; nd <= i__2; ++nd) {
4190 mma2cfv_(&jj, &minu, &maxu, nbpntv, &gssvtb[igsv *
4191 gssvtb_dim1], &chpair[(igsu + nd * chpair_dim2) *
4192 chpair_dim1], &chimpr[(igsu + nd * chimpr_dim2) *
4193 chimpr_dim1 + 1], &patjac[minu + (jj + nd *
4194 patjac_dim2) * patjac_dim1]);
4198 /* ----- Contribution of calculated terms to the approximation error */
4199 /* for terms (I,J) with MINU <= I <= MAXU, J fixe. */
4203 for (nd = 1; nd <= i__2; ++nd) {
4205 mma2er1_(ndjacu, ndjacv, &ndses, &minu, &maxu, &jj, &jj,
4206 iordru, iordrv, xmaxju, xmaxjv, &patjac[idim *
4207 patjac_dim2 * patjac_dim1], &vecerr[vecerr_dim1 + 1],
4208 &vecerr[nd + (vecerr_dim1 << 2)]);
4209 if (vecerr[nd + (vecerr_dim1 << 2)] > epsapr[nd]) {
4218 /* ******************************************************************
4220 /* ---------------------- Calculate the coeff of zone 2 -------------
4223 minu = (*iordru + 1) << 1;
4228 /* --> If zone 2 is empty, pass to zone 3. */
4229 /* VECERR(ND,2) was already set to zero. */
4234 /* ---------------- Calculate the terms connected to degree by U ------------
4238 for (nd = 1; nd <= i__1; ++nd) {
4240 for (kk = minu; kk <= i__2; ++kk) {
4242 mma2cfu_(&kk, nbpntu, nbpntv, &sosotb[nd * sosotb_dim2 *
4243 sosotb_dim1], &disotb[(nd * disotb_dim2 + 1) *
4244 disotb_dim1 + 1], &soditb[(nd * soditb_dim2 + 1) *
4245 soditb_dim1 + 1], &diditb[nd * diditb_dim2 *
4246 diditb_dim1], &gssutb[igsu * gssutb_dim1], &chpair[(
4247 igsu + nd * chpair_dim2) * chpair_dim1], &chimpr[(
4248 igsu + nd * chimpr_dim2) * chimpr_dim1 + 1]);
4254 /* ------------------- Calculate the coefficients of PATJAC ------------
4257 igsu = minu - i2rdu;
4259 for (jj = minv; jj <= i__1; ++jj) {
4262 for (nd = 1; nd <= i__2; ++nd) {
4263 mma2cfv_(&jj, &minu, &maxu, nbpntv, &gssvtb[igsv *
4264 gssvtb_dim1], &chpair[(igsu + nd * chpair_dim2) *
4265 chpair_dim1], &chimpr[(igsu + nd * chimpr_dim2) *
4266 chimpr_dim1 + 1], &patjac[minu + (jj + nd *
4267 patjac_dim2) * patjac_dim1]);
4273 /* -----Contribution of calculated terms to the approximation error */
4274 /* for terms (I,J) with MINU <= I <= MAXU, MINV <= J <= MAXV */
4278 for (nd = 1; nd <= i__1; ++nd) {
4280 mma2er1_(ndjacu, ndjacv, &ndses, &minu, &maxu, &minv, &maxv,
4281 iordru, iordrv, xmaxju, xmaxjv, &patjac[idim *
4282 patjac_dim2 * patjac_dim1], &vecerr[vecerr_dim1 + 1], &
4283 vecerr[nd + (vecerr_dim1 << 1)]);
4288 /* ******************************************************************
4290 /* ---------------------- Calculation of coeff of zone 3 -------------
4296 minv = (*iordrv + 1) << 1;
4299 /* -> If zone 3 is empty, pass to the test of cutting. */
4300 /* VECERR(ND,3) was already set to zero */
4305 /* ----------- The terms connected to the degree by U are already calculated -----
4307 /* ------------------- Calculation of coefficients of PATJAC ------------
4310 igsu = minu - i2rdu;
4312 for (jj = minv; jj <= i__1; ++jj) {
4315 for (nd = 1; nd <= i__2; ++nd) {
4316 mma2cfv_(&jj, &minu, &maxu, nbpntv, &gssvtb[igsv *
4317 gssvtb_dim1], &chpair[(igsu + nd * chpair_dim2) *
4318 chpair_dim1], &chimpr[(igsu + nd * chimpr_dim2) *
4319 chimpr_dim1 + 1], &patjac[minu + (jj + nd *
4320 patjac_dim2) * patjac_dim1]);
4326 /* ----- Contribution of calculated terms to the approximation error
4327 /* for terms (I,J) with MINU <= I <= MAXU, MINV <= J <= MAXV. */
4331 for (nd = 1; nd <= i__1; ++nd) {
4333 mma2er1_(ndjacu, ndjacv, &ndses, &minu, &maxu, &minv, &maxv,
4334 iordru, iordrv, xmaxju, xmaxjv, &patjac[idim *
4335 patjac_dim2 * patjac_dim1], &vecerr[vecerr_dim1 + 1], &
4336 vecerr[nd + vecerr_dim1 * 3]);
4341 /* ******************************************************************
4343 /* --------------------------- Tests of cutting ---------------------
4348 for (nd = 1; nd <= i__1; ++nd) {
4349 vaux[0] = vecerr[nd + (vecerr_dim1 << 1)];
4350 vaux[1] = vecerr[nd + (vecerr_dim1 << 2)];
4351 vaux[2] = vecerr[nd + vecerr_dim1 * 3];
4353 errmax[nd] = AdvApp2Var_MathBase::mzsnorm_(&ii, vaux);
4354 if (errmax[nd] > epsapr[nd]) {
4356 zv = AdvApp2Var_MathBase::mzsnorm_(&ii, vaux);
4357 zu = AdvApp2Var_MathBase::mzsnorm_(&ii, &vaux[1]);
4358 if (zu > epsapr[nd] && zv > epsapr[nd]) {
4370 /* ******************************************************************
4372 /* --- OK, the square is valid, the coeff of zone 1 are calculated
4375 minu = (*iordru + 1) << 1;
4377 minv = (*iordrv + 1) << 1;
4380 /* --> If zone 1 is empty, pass to the calculation of Max and Average error. */
4381 if (minu > maxu || minv > maxv) {
4385 /* ----------- The terms connected to degree by U are already calculated -----
4387 /* ------------------- Calculate the coefficients of PATJAC ------------
4390 igsu = minu - i2rdu;
4392 for (jj = minv; jj <= i__1; ++jj) {
4395 for (nd = 1; nd <= i__2; ++nd) {
4396 mma2cfv_(&jj, &minu, &maxu, nbpntv, &gssvtb[igsv *
4397 gssvtb_dim1], &chpair[(igsu + nd * chpair_dim2) *
4398 chpair_dim1], &chimpr[(igsu + nd * chimpr_dim2) *
4399 chimpr_dim1 + 1], &patjac[minu + (jj + nd *
4400 patjac_dim2) * patjac_dim1]);
4406 /* --------------- Now the degree is maximally lowered --------
4411 i__1 = 1, i__2 = (*iordru << 1) + 1, i__1 = advapp_max(i__1,i__2);
4412 minu = advapp_max(i__1,*ndminu);
4415 i__1 = 1, i__2 = (*iordrv << 1) + 1, i__1 = advapp_max(i__1,i__2);
4416 minv = advapp_max(i__1,*ndminv);
4420 for (nd = 1; nd <= i__1; ++nd) {
4422 if (maxu >= (*iordru + 1) << 1 && maxv >= (*iordrv + 1) << 1) {
4423 mma2er2_(ndjacu, ndjacv, &ndses, &minu, &maxu, &minv, &maxv,
4424 iordru, iordrv, xmaxju, xmaxjv, &patjac[idim *
4425 patjac_dim2 * patjac_dim1], &epsapr[nd], &vecerr[
4426 vecerr_dim1 + 1], &errmax[nd], &nu, &nv);
4434 /* --> Calculate the average error. */
4435 mma2moy_(ndjacu, ndjacv, &ndses, &nu1, ndjacu, &nv1, ndjacv,
4436 iordru, iordrv, &patjac[idim * patjac_dim2 * patjac_dim1],
4439 /* --> Set to 0.D0 the rejected coeffs. */
4440 i__2 = idim + ndses - 1;
4441 for (ii = idim; ii <= i__2; ++ii) {
4443 for (jj = nv1; jj <= i__3; ++jj) {
4445 for (kk = nu1; kk <= i__4; ++kk) {
4446 patjac[kk + (jj + ii * patjac_dim2) * patjac_dim1] =
4455 /* --> Return the nb of coeffs of approximation. */
4456 *ndegpu = advapp_max(*ndegpu,nu);
4457 *ndegpv = advapp_max(*ndegpv,nv);
4462 /* ******************************************************************
4464 /* -------------------- IT IS NOT POSSIBLE TO CUT -------------------
4466 /* ******************************************************************
4470 minu = (*iordru + 1) << 1;
4472 minv = (*iordrv + 1) << 1;
4475 /* ---------------- Calculate the terms connected to the degree by U ------------
4479 for (nd = 1; nd <= i__1; ++nd) {
4481 for (kk = minu; kk <= i__2; ++kk) {
4483 mma2cfu_(&kk, nbpntu, nbpntv, &sosotb[nd * sosotb_dim2 *
4484 sosotb_dim1], &disotb[(nd * disotb_dim2 + 1) *
4485 disotb_dim1 + 1], &soditb[(nd * soditb_dim2 + 1) *
4486 soditb_dim1 + 1], &diditb[nd * diditb_dim2 *
4487 diditb_dim1], &gssutb[igsu * gssutb_dim1], &chpair[(
4488 igsu + nd * chpair_dim2) * chpair_dim1], &chimpr[(
4489 igsu + nd * chimpr_dim2) * chimpr_dim1 + 1]);
4493 /* ---------------------- Calculate all coefficients -------
4496 igsu = minu - i2rdu;
4498 for (jj = minv; jj <= i__2; ++jj) {
4500 mma2cfv_(&jj, &minu, &maxu, nbpntv, &gssvtb[igsv *
4501 gssvtb_dim1], &chpair[(igsu + nd * chpair_dim2) *
4502 chpair_dim1], &chimpr[(igsu + nd * chimpr_dim2) *
4503 chimpr_dim1 + 1], &patjac[minu + (jj + nd *
4504 patjac_dim2) * patjac_dim1]);
4510 /* ----- Contribution of calculated terms to the approximation error
4511 /* for terms (I,J) with MINU <= I <= MAXU, MINV <= J <= MAXV */
4515 for (nd = 1; nd <= i__1; ++nd) {
4517 minu = (*iordru + 1) << 1;
4521 mma2er1_(ndjacu, ndjacv, &ndses, &minu, &maxu, &minv, &maxv,
4522 iordru, iordrv, xmaxju, xmaxjv, &patjac[idim *
4523 patjac_dim2 * patjac_dim1], &vecerr[vecerr_dim1 + 1], &
4527 minv = (*iordrv + 1) << 1;
4530 mma2er1_(ndjacu, ndjacv, &ndses, &minu, &maxu, &minv, &maxv,
4531 iordru, iordrv, xmaxju, xmaxjv, &patjac[idim *
4532 patjac_dim2 * patjac_dim1], &vecerr[vecerr_dim1 + 1],
4536 /* ---------------------------- IF ERRMAX > EPSAPR, stop --------
4539 if (errmax[nd] > epsapr[nd]) {
4544 /* ------------- Otherwise, try to remove again the coeff
4549 i__2 = 1, i__3 = (*iordru << 1) + 1, i__2 = advapp_max(i__2,i__3);
4550 minu = advapp_max(i__2,*ndminu);
4553 i__2 = 1, i__3 = (*iordrv << 1) + 1, i__2 = advapp_max(i__2,i__3);
4554 minv = advapp_max(i__2,*ndminv);
4556 if (maxu >= (*iordru + 1) << 1 && maxv >= (*iordrv + 1) << 1) {
4557 mma2er2_(ndjacu, ndjacv, &ndses, &minu, &maxu, &minv, &
4558 maxv, iordru, iordrv, xmaxju, xmaxjv, &patjac[
4559 idim * patjac_dim2 * patjac_dim1], &epsapr[nd], &
4560 vecerr[vecerr_dim1 + 1], &errmax[nd], &nu, &nv);
4567 /* --------------------- Calculate the average error -------------
4572 mma2moy_(ndjacu, ndjacv, &ndses, &nu1, ndjacu, &nv1, ndjacv,
4573 iordru, iordrv, &patjac[idim * patjac_dim2 * patjac_dim1],
4576 /* --------------------- Set to 0.D0 the rejected coeffs ----------
4579 i__2 = idim + ndses - 1;
4580 for (ii = idim; ii <= i__2; ++ii) {
4582 for (jj = nv1; jj <= i__3; ++jj) {
4584 for (kk = nu1; kk <= i__4; ++kk) {
4585 patjac[kk + (jj + ii * patjac_dim2) * patjac_dim1] =
4594 /* --------------- Return the nb of coeff of approximation ---
4597 *ndegpu = advapp_max(*ndegpu,nu);
4598 *ndegpv = advapp_max(*ndegpv,nv);
4606 /* ------------------------------ The end -------------------------------
4608 /* --> Error in inputs */
4613 /* --------- Management of cuts, it is required 0 < NUMDEC <= 5 -------
4616 /* --> Here it is possible and necessary to cut, choose by U if it is possible */
4618 if (*numdec <= 0 || *numdec > 5) {
4627 /* --> Here it is possible and necessary to cut, choose by U if it is possible */
4629 if (*numdec <= 0 || *numdec > 5) {
4638 /* --> Here it is possible and necessary to cut, choose by 4 if it is possible */
4640 if (*numdec <= 0 || *numdec > 5) {
4645 } else if (*numdec == 2 || *numdec == 4) {
4647 } else if (*numdec == 1 || *numdec == 3) {
4655 AdvApp2Var_SysBase::maermsg_("MMA2CE2", iercod, 7L);
4657 AdvApp2Var_SysBase::mgsomsg_("MMA2CE2", 7L);
4662 //=======================================================================
4663 //function : mma2cfu_
4665 //=======================================================================
4666 int mma2cfu_(integer *ndujac,
4678 /* System generated locals */
4679 integer sosotb_dim1, disotb_dim1, disotb_offset, soditb_dim1,
4680 soditb_offset, diditb_dim1, i__1, i__2;
4682 /* Local variables */
4684 integer nptu2, nptv2, ii, jj;
4685 doublereal bid0, bid1, bid2;
4687 /* **********************************************************************
4692 /* Calculate the terms connected to degree NDUJAC by U of the polynomial approximation */
4693 /* of function F(u,v), starting from its discretisation
4694 /* on the roots of Legendre polynom of degree */
4695 /* NBPNTU by U and NBPNTV by V. */
4699 /* FONCTION,APPROXIMATION,COEFFICIENT,POLYNOME */
4701 /* INPUT ARGUMENTSE : */
4702 /* ------------------ */
4703 /* NDUJAC: Fixed degree by U for which the terms */
4704 /* allowing to obtain the Legendre or Jacobi coeff*/
4705 /* of even or uneven degree by V are calculated. */
4706 /* NBPNTU: Degree of Legendre polynom on the roots which of */
4707 /* the coefficients of integration by U are calculated */
4708 /* by Gauss method. It is required that NBPNTU = 30, 40, 50 or 61. */
4709 /* NBPNTV: Degree of Legendre polynom on the roots which of */
4710 /* the coefficients of integration by V are calculated */
4711 /* by Gauss method. It is required that NBPNTV = 30, 40, 50 or 61. */
4712 /* SOSOTB: Table of F(ui,vj) + F(ui,-vj) + F(-ui,vj) + F(-ui,-vj) */
4713 /* with ui and vj positive roots of Legendre polynom */
4714 /* of degree NBPNTU and NBPNTV respectively. Moreover, */
4715 /* table SOSOTB(0,j) contains F(0,vj) + F(0,-vj), */
4716 /* table SOSOTB(i,0) contains F(ui,0) + F(-ui,0) and */
4717 /* SOSOTB(0,0) contains F(0,0). */
4718 /* DISOTB: Table of F(ui,vj) + F(ui,-vj) - F(-ui,vj) - F(-ui,-vj) */
4719 /* with ui and vj positive roots of Legendre polynom */
4720 /* of degree NBPNTU and NBPNTV respectively. */
4721 /* SODITB: Table of F(ui,vj) - F(ui,-vj) + F(-ui,vj) - F(-ui,-vj) */
4722 /* with ui and vj positive roots of Legendre polynom */
4723 /* of degree NBPNTU and NBPNTV respectively. */
4724 /* DIDITB: Table of F(ui,vj) - F(ui,-vj) - F(-ui,vj) + F(-ui,-vj) */
4725 /* avec ui and vj positive roots of Legendre polynom */
4726 /* of degree NBPNTU and NBPNTV respectively. Moreover, */
4727 /* table DIDITB(0,j) contains F(0,vj) - F(0,-vj), */
4728 /* and table DIDITB(i,0) contains F(ui,0) - F(-ui,0). */
4729 /* GSSUTB: Table of coefficients of integration by Gauss method */
4730 /* Gauss by U for fixed NDUJAC : i varies from 0 to NBPNTU/2. */
4732 /* OUTPUT ARGUMENTS : */
4733 /* ------------------- */
4734 /* CHPAIR: Table of terms connected to degree NDUJAC by U to calculate the */
4735 /* coeff. of the approximation of EVEN degree by V. */
4736 /* CHIMPR: Table of terms connected to degree NDUJAC by U to calculate */
4737 /* the coeff. of approximation of UNEVEN degree by V. */
4739 /* COMMONS USED : */
4740 /* ---------------- */
4742 /* REFERENCES CALLED : */
4743 /* ----------------------- */
4745 /* DESCRIPTION/NOTES/LIMITATIONS : */
4746 /* ----------------------------------- */
4750 /* **********************************************************************
4752 /* Name of the routine */
4755 /* --------------------------- Initialisations --------------------------
4758 /* Parameter adjustments */
4760 diditb_dim1 = *nbpntu / 2 + 1;
4761 soditb_dim1 = *nbpntu / 2;
4762 soditb_offset = soditb_dim1 + 1;
4763 soditb -= soditb_offset;
4764 disotb_dim1 = *nbpntu / 2;
4765 disotb_offset = disotb_dim1 + 1;
4766 disotb -= disotb_offset;
4767 sosotb_dim1 = *nbpntu / 2 + 1;
4770 ldbg = AdvApp2Var_SysBase::mnfndeb_() >= 3;
4772 AdvApp2Var_SysBase::mgenmsg_("MMA2CFU", 7L);
4775 nptu2 = *nbpntu / 2;
4776 nptv2 = *nbpntv / 2;
4778 /* **********************************************************************
4780 /* CALCULATE COEFFICIENTS BY U */
4782 /* ----------------- Calculate coefficients of even degree --------------
4785 if (*ndujac % 2 == 0) {
4787 for (jj = 1; jj <= i__1; ++jj) {
4791 for (ii = 1; ii <= i__2; ++ii) {
4793 bid1 += sosotb[ii + jj * sosotb_dim1] * bid0;
4794 bid2 += soditb[ii + jj * soditb_dim1] * bid0;
4802 /* --------------- Calculate coefficients of uneven degree ----------
4807 for (jj = 1; jj <= i__1; ++jj) {
4811 for (ii = 1; ii <= i__2; ++ii) {
4813 bid1 += disotb[ii + jj * disotb_dim1] * bid0;
4814 bid2 += diditb[ii + jj * diditb_dim1] * bid0;
4823 /* ------- Add terms connected to the supplementary root (0.D0) ------
4824 /* ----------- of Legendre polynom of uneven degree NBPNTU -----------
4826 /* --> Only even NDUJAC terms are modified as GSSUTB(0) = 0 */
4827 /* when NDUJAC is uneven. */
4829 if (*nbpntu % 2 != 0 && *ndujac % 2 == 0) {
4832 for (jj = 1; jj <= i__1; ++jj) {
4833 chpair[jj] += sosotb[jj * sosotb_dim1] * bid0;
4834 chimpr[jj] += diditb[jj * diditb_dim1] * bid0;
4839 /* ------ Calculate the terms connected to supplementary roots (0.D0) ------
4841 /* ----------- of Legendre polynom of uneven degree NBPNTV -----------
4844 if (*nbpntv % 2 != 0) {
4845 /* --> Only CHPAIR terms are calculated as GSSVTB(0,IH-IDEBV)=0
4847 /* when IH is uneven (see MMA2CFV). */
4849 if (*ndujac % 2 == 0) {
4852 for (ii = 1; ii <= i__1; ++ii) {
4853 bid1 += sosotb[ii] * gssutb[ii];
4860 for (ii = 1; ii <= i__1; ++ii) {
4861 bid1 += diditb[ii] * gssutb[ii];
4866 if (*nbpntu % 2 != 0) {
4867 chpair[0] += sosotb[0] * gssutb[0];
4871 /* ------------------------------ The end -------------------------------
4875 AdvApp2Var_SysBase::mgsomsg_("MMA2CFU", 7L);
4880 //=======================================================================
4881 //function : mma2cfv_
4883 //=======================================================================
4884 int mma2cfv_(integer *ndvjac,
4894 /* System generated locals */
4895 integer chpair_dim1, chpair_offset, chimpr_dim1, chimpr_offset,
4896 patjac_offset, i__1, i__2;
4898 /* Local variables */
4900 integer nptv2, ii, jj;
4903 /* **********************************************************************
4908 /* Calculate the coefficients of polynomial approximation of F(u,v)
4909 /* of degree NDVJAC by V and of degree by U varying from MINDGU to MAXDGU.
4914 /* FONCTION,APPROXIMATION,COEFFICIENT,POLYNOME */
4916 /* INPUT ARGUMENTS : */
4917 /* ------------------ */
4919 /* NDVJAC: Degree of the polynom of approximation by V. */
4920 /* The representation in the orthogonal base starts from degre 0.
4921 /* The polynomial base is the base of Jacobi of order -1 */
4922 /* (Legendre), 0, 1 or 2 */
4923 /* MINDGU: Degree minimum by U of coeff. to calculate. */
4924 /* MAXDGU: Degree maximum by U of coeff. to calculate. */
4925 /* NBPNTV: Degree of the Legendre polynom on the roots which of */
4926 /* the coefficients of integration by V are calculated */
4927 /* by Gauss method. It is reqired that NBPNTV = 30, 40, 50 or 61 and NDVJAC < NBPNTV. */
4928 /* GSSVTB: Table of coefficients of integration by Gauss method */
4929 /* by V for NDVJAC fixed: j varies from 0 to NBPNTV/2. */
4930 /* CHPAIR: Table of terms connected to degrees from MINDGU to MAXDGU by U to
4931 /* calculate the coeff. of approximation of EVEN degree NDVJAC by V. */
4932 /* CHIMPR: Table of terms connected to degrees from MINDGU to MAXDGU by U to
4933 /* calculate the coeff. of approximation of UNEVEN degree NDVJAC by V. */
4935 /* OUTPUT ARGUMENTS : */
4936 /* ------------------- */
4937 /* PATJAC: Table of coefficients by U of the polynom of approximation */
4938 /* P(u,v) of degree MINDGU to MAXDGU by U and NDVJAC by V. */
4940 /* COMMONS USED : */
4941 /* -------------- */
4943 /* REFERENCES CALLED : */
4944 /* --------------------- */
4946 /* DESCRIPTION/NOTES/LIMITATIONS : */
4947 /* ------------------------------- */
4949 /* **********************************************************************
4951 /* Name of the routine */
4954 /* --------------------------- Initialisations --------------------------
4957 /* Parameter adjustments */
4958 patjac_offset = *mindgu;
4959 patjac -= patjac_offset;
4960 chimpr_dim1 = *nbpntv / 2;
4961 chimpr_offset = chimpr_dim1 * *mindgu + 1;
4962 chimpr -= chimpr_offset;
4963 chpair_dim1 = *nbpntv / 2 + 1;
4964 chpair_offset = chpair_dim1 * *mindgu;
4965 chpair -= chpair_offset;
4968 ldbg = AdvApp2Var_SysBase::mnfndeb_() >= 3;
4970 AdvApp2Var_SysBase::mgenmsg_("MMA2CFV", 7L);
4972 nptv2 = *nbpntv / 2;
4974 /* --------- Calculate the coefficients for even degree NDVJAC ----------
4977 if (*ndvjac % 2 == 0) {
4979 for (ii = *mindgu; ii <= i__1; ++ii) {
4982 for (jj = 1; jj <= i__2; ++jj) {
4983 bid1 += chpair[jj + ii * chpair_dim1] * gssvtb[jj];
4990 /* -------- Calculate the coefficients for uneven degree NDVJAC -----
4995 for (ii = *mindgu; ii <= i__1; ++ii) {
4998 for (jj = 1; jj <= i__2; ++jj) {
4999 bid1 += chimpr[jj + ii * chimpr_dim1] * gssvtb[jj];
5007 /* ------- Add terms connected to the supplementary root (0.D0) ----- */
5008 /* --------of the Legendre polynom of uneven degree NBPNTV --------- */
5010 if (*nbpntv % 2 != 0 && *ndvjac % 2 == 0) {
5013 for (ii = *mindgu; ii <= i__1; ++ii) {
5014 patjac[ii] += bid1 * chpair[ii * chpair_dim1];
5019 /* ------------------------------ The end -------------------------------
5023 AdvApp2Var_SysBase::mgsomsg_("MMA2CFV", 7L);
5028 //=======================================================================
5029 //function : mma2ds1_
5031 //=======================================================================
5032 int AdvApp2Var_ApproxF2var::mma2ds1_(integer *ndimen,
5035 const AdvApp2Var_EvaluatorFunc2Var& foncnp,
5050 /* System generated locals */
5051 integer sosotb_dim1, sosotb_dim2, sosotb_offset, disotb_dim1, disotb_dim2,
5052 disotb_offset, soditb_dim1, soditb_dim2, soditb_offset,
5053 diditb_dim1, diditb_dim2, diditb_offset, fpntab_dim1,
5054 fpntab_offset, i__1;
5056 /* Local variables */
5058 integer ibid1, ibid2, iuouv, nd;
5061 /* **********************************************************************
5066 /* Discretisation of function F(u,v) on the roots of Legendre polynoms. */
5070 /* FONCTION&,DISCRETISATION,&POINT */
5072 /* INPUT ARGUMENTS : */
5073 /* ------------------ */
5074 /* NDIMEN: Dimension of the space. */
5075 /* UINTFN: Limits of the interval of definition by u of the function */
5076 /* to be processed: (UINTFN(1),UINTFN(2)). */
5077 /* VINTFN: Limits of the interval of definition by v of the function */
5078 /* to be processed: (VINTFN(1),VINTFN(2)). */
5079 /* FONCNP: The NAME of the non-polynomial function to be processed. */
5080 /* NBPNTU: The degree of Legendre polynom on the roots which of */
5081 /* FONCNP is discretized by u. */
5082 /* NBPNTV: The degree of Legendre polynom on the roots which of */
5083 /* FONCNP is discretized by v. */
5084 /* UROOTB: Table of STRICTLY POSITIVE roots of the polynom */
5085 /* of Legendre of degree NBPNTU defined on (-1,1). */
5086 /* VROOTB: Table of STRICTLY POSITIVE roots of the polynom */
5087 /* of Legendre of degree NBPNTV defined on (-1,1). */
5088 /* ISOFAV: Shows the type of iso of F(u,v) to be extracted to improve */
5089 /* the rapidity of calculation (has no influence on the form */
5091 /* = 1, shows that it is necessary to calculate the points of F(u,v) */
5092 /* with fixed u (with NBPNTV values different from v). */
5093 /* = 2, shows that it is necessaty to calculate the points of F(u,v) */
5094 /* with fixed v (with NBPNTU values different from u). */
5095 /* SOSOTB: Preinitialized table (input/output argument). */
5096 /* DISOTB: Preinitialized table (input/output argument). */
5097 /* SODITB: Preinitialized table (input/output argument). */
5098 /* DIDITB: Preinitialized table (input/output argument). */
5100 /* OUTPUT ARGUMENTS : */
5101 /* ------------------- */
5102 /* SOSOTB: Table where the terms */
5103 /* F(ui,vj) + F(ui,-vj) + F(-ui,vj) + F(-ui,-vj) */
5104 /* are added with ui and vj positive roots of Legendre polynom */
5105 /* of degree NBPNTU and NBPNTV respectively. */
5106 /* DISOTB: Table where the terms */
5107 /* F(ui,vj) + F(ui,-vj) - F(-ui,vj) - F(-ui,-vj) */
5108 /* are added with ui and vj positive roots of Legendre polynom */
5109 /* of degree NBPNTU and NBPNTV respectively. */
5110 /* SODITB: Table where the terms */
5111 /* F(ui,vj) - F(ui,-vj) + F(-ui,vj) - F(-ui,-vj) */
5112 /* are added with ui and vj positive roots of Legendre polynom */
5113 /* of degree NBPNTU and NBPNTV respectively. */
5114 /* DIDITB: Table where the terms */
5115 /* F(ui,vj) - F(ui,-vj) - F(-ui,vj) + F(-ui,-vj) */
5116 /* are added with ui and vj positive roots of Legendre polynom */
5117 /* of degree NBPNTU and NBPNTV respectively. */
5118 /* FPNTAB: Auxiliary table. */
5119 /* TTABLE: Auxiliary table. */
5120 /* IERCOD: Error code >100 Pb in the evaluation of FONCNP, */
5121 /* the returned error code is equal to error code of FONCNP + 100. */
5123 /* COMMONS USED : */
5124 /* ---------------- */
5126 /* REFERENCES CALLED : */
5127 /* --------------------- */
5129 /* DESCRIPTION/NOTES/LIMITATIONS : */
5130 /* ----------------------------------- */
5131 /* --> The external function created by the caller of MA2F1K, MA2FDK */
5132 /* where MA2FXK should be in the following form : */
5133 /* SUBROUTINE FONCNP(NDIMEN,UINTFN,VINTFN,ISOFAV,TCONST,NBPTAB */
5134 /* ,TTABLE,IDERIU,IDERIV,PPNTAB,IERCOD) */
5135 /* with the following input arguments : */
5136 /* - NDIMEN is integer defined as the sum of dimensions of */
5137 /* sub-spaces (i.e. total dimension of the problem). */
5138 /* - UINTFN(2) is a table of 2 reals containing the interval */
5139 /* by u where the function to be approximated is defined */
5140 /* (so it is equal to UIFONC). */
5141 /* - VINTFN(2) is a table of 2 reals containing the interval */
5142 /* by v where the function to be approximated is defined */
5143 /* (so it is equal to VIFONC). */
5144 /* - ISOFAV, is 1 if it is necessary to calculate points with constant u, */
5145 /* is 2 if it is necessary to calculate points with constant v. */
5146 /* Any other value is an error. */
5147 /* - TCONST, real, value of the fixed parameter. Takes values */
5148 /* in (UIFONC(1),UIFONC(2)) if ISOFAV = 1 or */
5149 /* ins (VIFONC(1),VIFONC(2)) if ISOFAV = 2. */
5150 /* - NBPTAB, integer. Shows the number of points to be calculated. */
5151 /* - TTABLE, a table of reals NBPTAB. These are the values of */
5152 /* 'free' parameter of discretization (v if IISOFAV=1, */
5153 /* u if IISOFAV=2). */
5154 /* - IDERIU, integer, takes values between 0 (position) */
5155 /* and IORDRE(1) (partial derivative of the function by u */
5156 /* of order IORDRE(1) if IORDRE(1) > 0). */
5157 /* - IDERIV, integer, takes values between 0 (position) */
5158 /* and IORDRE(2) (partial derivative of the function by v */
5159 /* of order IORDRE(2) if IORDRE(2) > 0). */
5160 /* If IDERIU=i and IDERIV=j, FONCNP should calculate the */
5161 /* points of the derivative : */
5168 /* and the output arguments aret : */
5169 /* - FPNTAB(NDIMEN,NBPTAB) contains, at output, the table of */
5170 /* NBPTAB points calculated in FONCNP. */
5171 /* - IERCOD is, at output the error code of FONCNP. This code */
5172 /* (integer) should be strictly positive if there is a problem. */
5174 /* The input arguments SHOULD NOT be modified under FONCNP.
5177 /* -->As FONCNP is not forcedly defined in (-1,1)*(-1,1), the */
5178 /* values of UROOTB and VROOTB are consequently modified. */
5180 /* -->The results of discretisation are ranked in 4 tables */
5181 /* SOSOTB, DISOTB, SODITB and DIDITB to earn time */
5182 /* during the calculation of coefficients of the polynom of approximation. */
5184 /* When NBPNTU is uneven : */
5185 /* table SOSOTB(0,j) contains F(0,vj) + F(0,-vj), */
5186 /* table DIDITB(0,j) contains F(0,vj) - F(0,-vj), */
5187 /* When NBPNTV is uneven : */
5188 /* table SOSOTB(i,0) contains F(ui,0) + F(-ui,0), */
5189 /* table DIDITB(i,0) contains F(ui,0) - F(-ui,0), */
5190 /* When NBPNTU and NBPNTV are uneven : */
5191 /* term SOSOTB(0,0) contains F(0,0). */
5194 /* **********************************************************************
5196 /* Name of the routine */
5199 /* --------------------------- Initialization --------------------------
5202 /* Parameter adjustments */
5203 fpntab_dim1 = *ndimen;
5204 fpntab_offset = fpntab_dim1 + 1;
5205 fpntab -= fpntab_offset;
5209 diditb_dim1 = *nbpntu / 2 + 1;
5210 diditb_dim2 = *nbpntv / 2 + 1;
5211 diditb_offset = diditb_dim1 * diditb_dim2;
5212 diditb -= diditb_offset;
5213 soditb_dim1 = *nbpntu / 2;
5214 soditb_dim2 = *nbpntv / 2;
5215 soditb_offset = soditb_dim1 * (soditb_dim2 + 1) + 1;
5216 soditb -= soditb_offset;
5217 disotb_dim1 = *nbpntu / 2;
5218 disotb_dim2 = *nbpntv / 2;
5219 disotb_offset = disotb_dim1 * (disotb_dim2 + 1) + 1;
5220 disotb -= disotb_offset;
5221 sosotb_dim1 = *nbpntu / 2 + 1;
5222 sosotb_dim2 = *nbpntv / 2 + 1;
5223 sosotb_offset = sosotb_dim1 * sosotb_dim2;
5224 sosotb -= sosotb_offset;
5229 ldbg = AdvApp2Var_SysBase::mnfndeb_() >= 3;
5231 AdvApp2Var_SysBase::mgenmsg_("MMA2DS1", 7L);
5234 if (*isofav < 1 || *isofav > 2) {
5240 /* **********************************************************************
5242 /* --------- Discretization by U on the roots of the polynom of ------ */
5243 /* --------------- Legendre of degree NBPNTU, iso-V by iso-V --------- */
5244 /* **********************************************************************
5248 mma2ds2_(ndimen, &uintfn[1], &vintfn[1], foncnp, nbpntu, nbpntv, &
5249 urootb[1], &vrootb[1], &iuouv, &sosotb[sosotb_offset], &
5250 disotb[disotb_offset], &soditb[soditb_offset], &diditb[
5251 diditb_offset], &fpntab[fpntab_offset], &ttable[1], iercod);
5253 /* ******************************************************************
5255 /* --------- Discretization by V on the roots of the polynom of ------ */
5256 /* --------------- Legendre of degree NBPNTV, iso-V by iso-V --------- */
5257 /* ******************************************************************
5261 /* --> Inversion of indices of tables */
5263 for (nd = 1; nd <= i__1; ++nd) {
5264 isz1 = *nbpntu / 2 + 1;
5265 isz2 = *nbpntv / 2 + 1;
5266 AdvApp2Var_MathBase::mmfmtb1_(&isz1, &sosotb[nd * sosotb_dim2 * sosotb_dim1], &isz1, &
5267 isz2, &isz2, &sosotb[nd * sosotb_dim2 * sosotb_dim1], &
5268 ibid1, &ibid2, iercod);
5272 AdvApp2Var_MathBase::mmfmtb1_(&isz1, &diditb[nd * diditb_dim2 * diditb_dim1], &isz1, &
5273 isz2, &isz2, &diditb[nd * diditb_dim2 * diditb_dim1], &
5274 ibid1, &ibid2, iercod);
5280 AdvApp2Var_MathBase::mmfmtb1_(&isz1, &soditb[(nd * soditb_dim2 + 1) * soditb_dim1 + 1],
5281 &isz1, &isz2, &isz2, &soditb[(nd * soditb_dim2 + 1) *
5282 soditb_dim1 + 1], &ibid1, &ibid2, iercod);
5286 AdvApp2Var_MathBase::mmfmtb1_(&isz1, &disotb[(nd * disotb_dim2 + 1) * disotb_dim1 + 1],
5287 &isz1, &isz2, &isz2, &disotb[(nd * disotb_dim2 + 1) *
5288 disotb_dim1 + 1], &ibid1, &ibid2, iercod);
5295 mma2ds2_(ndimen, &vintfn[1], &uintfn[1], foncnp, nbpntv, nbpntu, &
5296 vrootb[1], &urootb[1], &iuouv, &sosotb[sosotb_offset], &
5297 soditb[soditb_offset], &disotb[disotb_offset], &diditb[
5298 diditb_offset], &fpntab[fpntab_offset], &ttable[1], iercod);
5299 /* --> Inversion of indices of tables */
5301 for (nd = 1; nd <= i__1; ++nd) {
5302 isz1 = *nbpntv / 2 + 1;
5303 isz2 = *nbpntu / 2 + 1;
5304 AdvApp2Var_MathBase::mmfmtb1_(&isz1, &sosotb[nd * sosotb_dim2 * sosotb_dim1], &isz1, &
5305 isz2, &isz2, &sosotb[nd * sosotb_dim2 * sosotb_dim1], &
5306 ibid1, &ibid2, iercod);
5310 AdvApp2Var_MathBase::mmfmtb1_(&isz1, &diditb[nd * diditb_dim2 * diditb_dim1], &isz1, &
5311 isz2, &isz2, &diditb[nd * diditb_dim2 * diditb_dim1], &
5312 ibid1, &ibid2, iercod);
5318 AdvApp2Var_MathBase::mmfmtb1_(&isz1, &soditb[(nd * soditb_dim2 + 1) * soditb_dim1 + 1],
5319 &isz1, &isz2, &isz2, &soditb[(nd * soditb_dim2 + 1) *
5320 soditb_dim1 + 1], &ibid1, &ibid2, iercod);
5324 AdvApp2Var_MathBase::mmfmtb1_(&isz1, &disotb[(nd * disotb_dim2 + 1) * disotb_dim1 + 1],
5325 &isz1, &isz2, &isz2, &disotb[(nd * disotb_dim2 + 1) *
5326 disotb_dim1 + 1], &ibid1, &ibid2, iercod);
5334 /* ------------------------------ The end -------------------------------
5340 AdvApp2Var_SysBase::maermsg_("MMA2DS1", iercod, 7L);
5343 AdvApp2Var_SysBase::mgsomsg_("MMA2DS1", 7L);
5348 //=======================================================================
5349 //function : mma2ds2_
5351 //=======================================================================
5352 int mma2ds2_(integer *ndimen,
5355 const AdvApp2Var_EvaluatorFunc2Var& foncnp,
5371 /* System generated locals */
5372 integer sosotb_dim1, sosotb_dim2, sosotb_offset, disotb_dim1, disotb_dim2,
5373 disotb_offset, soditb_dim1, soditb_dim2, soditb_offset,
5374 diditb_dim1, diditb_dim2, diditb_offset, fpntab_dim1,
5375 fpntab_offset, i__1, i__2, i__3;
5377 /* Local variables */
5380 doublereal alinu, blinu, alinv, blinv, tcons;
5381 doublereal dbfn1[2], dbfn2[2];
5382 integer nuroo, nvroo, id, iu, iv;
5386 /* **********************************************************************
5391 /* Discretization of function F(u,v) on the roots of polynoms of Legendre. */
5395 /* FONCTION&,DISCRETISATION,&POINT */
5397 /* INPUT ARGUMENTS : */
5398 /* ------------------ */
5399 /* NDIMEN: Dimension of the space. */
5400 /* UINTFN: Limits of the interval of definition by u of the function */
5401 /* to be processed: (UINTFN(1),UINTFN(2)). */
5402 /* VINTFN: Limits of the interval of definition by v of the function */
5403 /* to be processed: (VINTFN(1),VINTFN(2)). */
5404 /* FONCNP: The NAME of the non-polynomial function to be processed. */
5405 /* NBPNTU: The degree of Legendre polynom on the roots which of */
5406 /* FONCNP is discretized by u. */
5407 /* NBPNTV: The degree of Legendre polynom on the roots which of */
5408 /* FONCNP is discretized by v. */
5409 /* UROOTB: Table of STRICTLY POSITIVE roots of the polynom */
5410 /* of Legendre of degree NBPNTU defined on (-1,1). */
5411 /* VROOTB: Table of STRICTLY POSITIVE roots of the polynom */
5412 /* of Legendre of degree NBPNTV defined on (-1,1). */
5413 /* IIUOUV: Shows the type of iso of F(u,v) tom be extracted to improve the */
5414 /* rapidity of calculation (has no influence on the form of result) */
5415 /* = 1, shows that it is necessary to calculate the points of F(u,v) */
5416 /* with fixed u (so with NBPNTV values different from v). */
5417 /* = 2, shows that it is necessary to calculate the points of F(u,v) */
5418 /* with fixed v (so with NBPNTV values different from u). */
5419 /* SOSOTB: Preinitialized table (input/output argument). */
5420 /* DISOTB: Preinitialized table (input/output argument). */
5421 /* SODITB: Preinitialized table (input/output argument). */
5422 /* DIDITB: Preinitialized table (input/output argument). */
5424 /* OUTPUT ARGUMENTS : */
5425 /* ------------------- */
5426 /* SOSOTB: Table where the terms */
5427 /* F(ui,vj) + F(ui,-vj) + F(-ui,vj) + F(-ui,-vj) */
5428 /* are added with ui and vj positive roots of Legendre polynom */
5429 /* of degree NBPNTU and NBPNTV respectively. */
5430 /* DISOTB: Table where the terms */
5431 /* F(ui,vj) + F(ui,-vj) - F(-ui,vj) - F(-ui,-vj) */
5432 /* are added with ui and vj positive roots of Legendre polynom */
5433 /* of degree NBPNTU and NBPNTV respectively. */
5434 /* SODITB: Table where the terms */
5435 /* F(ui,vj) - F(ui,-vj) + F(-ui,vj) - F(-ui,-vj) */
5436 /* are added with ui and vj positive roots of Legendre polynom */
5437 /* of degree NBPNTU and NBPNTV respectively. */
5438 /* DIDITB: Table where the terms */
5439 /* F(ui,vj) - F(ui,-vj) - F(-ui,vj) + F(-ui,-vj) */
5440 /* are added with ui and vj positive roots of Legendre polynom */
5441 /* of degree NBPNTU and NBPNTV respectively. */
5442 /* FPNTAB: Auxiliary table. */
5443 /* TTABLE: Auxiliary table. */
5444 /* IERCOD: Error code >100 Pb in the evaluation of FONCNP, */
5445 /* the returned error code is equal to error code of FONCNP + 100. */
5447 /* COMMONS USED : */
5448 /* ---------------- */
5450 /* REFERENCES CALLED : */
5451 /* --------------------- */
5453 /* DESCRIPTION/NOTES/LIMITATIONS : */
5454 /* ----------------------------------- */
5455 /* --> The external function created by the caller of MA2F1K, MA2FDK */
5456 /* where MA2FXK should be in the following form : */
5457 /* SUBROUTINE FONCNP(NDIMEN,UINTFN,VINTFN,IIIUOUV,TCONST,NBPTAB */
5458 /* ,TTABLE,IDERIU,IDERIV,PPNTAB,IERCOD) */
5459 /* with the following input arguments : */
5460 /* - NDIMEN is integer defined as the sum of dimensions of */
5461 /* sub-spaces (i.e. total dimension of the problem). */
5462 /* - UINTFN(2) is a table of 2 reals containing the interval */
5463 /* by u where the function to be approximated is defined */
5464 /* (so it is equal to UIFONC). */
5465 /* - VINTFN(2) is a table of 2 reals containing the interval */
5466 /* by v where the function to be approximated is defined */
5467 /* (so it is equal to VIFONC). */
5468 /* - IIIUOUV, is 1 if it is necessary to calculate points with constant u, */
5469 /* is 2 if it is necessary to calculate points with constant v. */
5470 /* Any other value is an error. */
5471 /* - TCONST, real, value of the fixed parameter. Takes values */
5472 /* in (UIFONC(1),UIFONC(2)) if ISOFAV = 1 or */
5473 /* ins (VIFONC(1),VIFONC(2)) if ISOFAV = 2. */
5474 /* - NBPTAB, integer. Shows the number of points to be calculated. */
5475 /* - TTABLE, a table of reals NBPTAB. These are the values of */
5476 /* 'free' parameter of discretization (v if IIIUOUV=1, */
5477 /* u if IIIUOUV=2). */
5478 /* - IDERIU, integer, takes values between 0 (position) */
5479 /* and IORDRE(1) (partial derivative of the function by u */
5480 /* of order IORDRE(1) if IORDRE(1) > 0). */
5481 /* - IDERIV, integer, takes values between 0 (position) */
5482 /* and IORDRE(2) (partial derivative of the function by v */
5483 /* of order IORDRE(2) if IORDRE(2) > 0). */
5484 /* If IDERIU=i and IDERIV=j, FONCNP should calculate the */
5485 /* points of the derivative : */
5492 /* and the output arguments aret : */
5493 /* - FPNTAB(NDIMEN,NBPTAB) contains, at output, the table of */
5494 /* NBPTAB points calculated in FONCNP. */
5495 /* - IERCOD is, at output the error code of FONCNP. This code */
5496 /* (integer) should be strictly positive if there is a problem. */
5498 /* The input arguments SHOULD NOT be modified under FONCNP.
5501 /* -->As FONCNP is not forcedly defined in (-1,1)*(-1,1), the */
5502 /* values of UROOTB and VROOTB are consequently modified. */
5504 /* -->The results of discretisation are ranked in 4 tables */
5505 /* SOSOTB, DISOTB, SODITB and DIDITB to earn time */
5506 /* during the calculation of coefficients of the polynom of approximation. */
5508 /* When NBPNTU is uneven : */
5509 /* table SOSOTB(0,j) contains F(0,vj) + F(0,-vj), */
5510 /* table DIDITB(0,j) contains F(0,vj) - F(0,-vj), */
5511 /* When NBPNTV is uneven : */
5512 /* table SOSOTB(i,0) contains F(ui,0) + F(-ui,0), */
5513 /* table DIDITB(i,0) contains F(ui,0) - F(-ui,0), */
5514 /* When NBPNTU and NBPNTV are uneven : */
5515 /* term SOSOTB(0,0) contains F(0,0). */
5517 /* ATTENTION: These 4 tables are filled by varying the */
5518 /* 1st index first. So, the discretizations */
5519 /* of F(...,t) (for IIUOUV = 2) or of F(t,...) (IIUOUV = 1) */
5520 /* are stored in SOSOTB(...,t), SODITB(...,t), etc... */
5521 /* (this allows to gain important time). */
5522 /* It is required that the caller, in case of IIUOUV=1, */
5523 /* invert the roles of u and v, of SODITB and DISOTB BEFORE the */
5526 /* **********************************************************************
5529 /* Name of the routine */
5531 /* --> Indices of loops. */
5533 /* --------------------------- Initialization --------------------------
5536 /* Parameter adjustments */
5540 fpntab_dim1 = *ndimen;
5541 fpntab_offset = fpntab_dim1 + 1;
5542 fpntab -= fpntab_offset;
5544 diditb_dim1 = *nbpntu / 2 + 1;
5545 diditb_dim2 = *nbpntv / 2 + 1;
5546 diditb_offset = diditb_dim1 * diditb_dim2;
5547 diditb -= diditb_offset;
5548 soditb_dim1 = *nbpntu / 2;
5549 soditb_dim2 = *nbpntv / 2;
5550 soditb_offset = soditb_dim1 * (soditb_dim2 + 1) + 1;
5551 soditb -= soditb_offset;
5552 disotb_dim1 = *nbpntu / 2;
5553 disotb_dim2 = *nbpntv / 2;
5554 disotb_offset = disotb_dim1 * (disotb_dim2 + 1) + 1;
5555 disotb -= disotb_offset;
5556 sosotb_dim1 = *nbpntu / 2 + 1;
5557 sosotb_dim2 = *nbpntv / 2 + 1;
5558 sosotb_offset = sosotb_dim1 * sosotb_dim2;
5559 sosotb -= sosotb_offset;
5563 ldbg = AdvApp2Var_SysBase::mnfndeb_() >= 3;
5565 AdvApp2Var_SysBase::mgenmsg_("MMA2DS2", 7L);
5569 alinu = (uintfn[2] - uintfn[1]) / 2.;
5570 blinu = (uintfn[2] + uintfn[1]) / 2.;
5571 alinv = (vintfn[2] - vintfn[1]) / 2.;
5572 blinv = (vintfn[2] + vintfn[1]) / 2.;
5575 dbfn1[0] = vintfn[1];
5576 dbfn1[1] = vintfn[2];
5577 dbfn2[0] = uintfn[1];
5578 dbfn2[1] = uintfn[2];
5580 dbfn1[0] = uintfn[1];
5581 dbfn1[1] = uintfn[2];
5582 dbfn2[0] = vintfn[1];
5583 dbfn2[1] = vintfn[2];
5586 /* **********************************************************************
5588 /* -------- Discretization by U on the roots of Legendre polynom -------- */
5589 /* ---------------- of degree NBPNTU, with Vj fixed -------------------- */
5590 /* **********************************************************************
5593 nuroo = *nbpntu / 2;
5594 nvroo = *nbpntv / 2;
5595 jdec = (*nbpntu + 1) / 2;
5597 /* ----------- Loading of parameters of discretization by U ------------- */
5600 for (iu = 1; iu <= i__1; ++iu) {
5601 ttable[iu] = blinu + alinu * urootb[iu];
5605 /* -------------- For Vj fixed, negative root of Legendre ------------- */
5608 for (iv = 1; iv <= i__1; ++iv) {
5609 tcons = blinv + alinv * vrootb[iv];
5610 (*const_cast <AdvApp2Var_EvaluatorFunc2Var*> (&foncnp)).Evaluate (
5611 ndimen, dbfn1, dbfn2, iiuouv, &tcons, nbpntu, &
5612 ttable[1], &c__0, &c__0, &fpntab[fpntab_offset], iercod);
5617 for (id = 1; id <= i__2; ++id) {
5619 for (iu = 1; iu <= i__3; ++iu) {
5620 up = fpntab[id + (iu + jdec) * fpntab_dim1];
5621 um = fpntab[id + (nuroo - iu + 1) * fpntab_dim1];
5622 sosotb[iu + (nvroo - iv + 1 + id * sosotb_dim2) * sosotb_dim1]
5623 = sosotb[iu + (nvroo - iv + 1 + id * sosotb_dim2) *
5624 sosotb_dim1] + up + um;
5625 disotb[iu + (nvroo - iv + 1 + id * disotb_dim2) * disotb_dim1]
5626 = disotb[iu + (nvroo - iv + 1 + id * disotb_dim2) *
5627 disotb_dim1] + up - um;
5628 soditb[iu + (nvroo - iv + 1 + id * soditb_dim2) * soditb_dim1]
5629 = soditb[iu + (nvroo - iv + 1 + id * soditb_dim2) *
5630 soditb_dim1] - up - um;
5631 diditb[iu + (nvroo - iv + 1 + id * diditb_dim2) * diditb_dim1]
5632 = diditb[iu + (nvroo - iv + 1 + id * diditb_dim2) *
5633 diditb_dim1] - up + um;
5636 if (*nbpntu % 2 != 0) {
5637 up = fpntab[id + jdec * fpntab_dim1];
5638 sosotb[(nvroo - iv + 1 + id * sosotb_dim2) * sosotb_dim1] +=
5640 diditb[(nvroo - iv + 1 + id * diditb_dim2) * diditb_dim1] -=
5648 /* --------- For Vj = 0 (uneven NBPNTV), discretization by U ----------- */
5650 if (*nbpntv % 2 != 0) {
5652 (*const_cast <AdvApp2Var_EvaluatorFunc2Var*> (&foncnp)).Evaluate (
5653 ndimen, dbfn1, dbfn2, iiuouv, &tcons, nbpntu, &
5654 ttable[1], &c__0, &c__0, &fpntab[fpntab_offset], iercod);
5659 for (id = 1; id <= i__1; ++id) {
5661 for (iu = 1; iu <= i__2; ++iu) {
5662 up = fpntab[id + (jdec + iu) * fpntab_dim1];
5663 um = fpntab[id + (nuroo - iu + 1) * fpntab_dim1];
5664 sosotb[iu + id * sosotb_dim2 * sosotb_dim1] = sosotb[iu + id *
5665 sosotb_dim2 * sosotb_dim1] + up + um;
5666 diditb[iu + id * diditb_dim2 * diditb_dim1] = diditb[iu + id *
5667 diditb_dim2 * diditb_dim1] + up - um;
5670 if (*nbpntu % 2 != 0) {
5671 up = fpntab[id + jdec * fpntab_dim1];
5672 sosotb[id * sosotb_dim2 * sosotb_dim1] += up;
5678 /* -------------- For Vj fixed, positive root of Legendre ------------- */
5681 for (iv = 1; iv <= i__1; ++iv) {
5682 tcons = alinv * vrootb[(*nbpntv + 1) / 2 + iv] + blinv;
5683 (*const_cast <AdvApp2Var_EvaluatorFunc2Var*> (&foncnp)).Evaluate (
5684 ndimen, dbfn1, dbfn2, iiuouv, &tcons, nbpntu, &
5685 ttable[1], &c__0, &c__0, &fpntab[fpntab_offset], iercod);
5690 for (id = 1; id <= i__2; ++id) {
5692 for (iu = 1; iu <= i__3; ++iu) {
5693 up = fpntab[id + (iu + jdec) * fpntab_dim1];
5694 um = fpntab[id + (nuroo - iu + 1) * fpntab_dim1];
5695 sosotb[iu + (iv + id * sosotb_dim2) * sosotb_dim1] = sosotb[
5696 iu + (iv + id * sosotb_dim2) * sosotb_dim1] + up + um;
5697 disotb[iu + (iv + id * disotb_dim2) * disotb_dim1] = disotb[
5698 iu + (iv + id * disotb_dim2) * disotb_dim1] + up - um;
5699 soditb[iu + (iv + id * soditb_dim2) * soditb_dim1] = soditb[
5700 iu + (iv + id * soditb_dim2) * soditb_dim1] + up + um;
5701 diditb[iu + (iv + id * diditb_dim2) * diditb_dim1] = diditb[
5702 iu + (iv + id * diditb_dim2) * diditb_dim1] + up - um;
5705 if (*nbpntu % 2 != 0) {
5706 up = fpntab[id + jdec * fpntab_dim1];
5707 sosotb[(iv + id * sosotb_dim2) * sosotb_dim1] += up;
5708 diditb[(iv + id * diditb_dim2) * diditb_dim1] += up;
5715 /* ------------------------------ The end -------------------------------
5721 AdvApp2Var_SysBase::maermsg_("MMA2DS2", iercod, 7L);
5724 AdvApp2Var_SysBase::mgsomsg_("MMA2DS2", 7L);
5729 //=======================================================================
5730 //function : mma2er1_
5732 //=======================================================================
5733 int mma2er1_(integer *ndjacu,
5749 /* System generated locals */
5750 integer patjac_dim1, patjac_dim2, patjac_offset, i__1, i__2, i__3;
5753 /* Local variables */
5758 doublereal bid0, bid1;
5760 /* **********************************************************************
5765 /* Calculate max approximation error done when */
5766 /* the coefficients of PATJAC such that the degree by U varies between */
5767 /* MINDGU and MAXDGU and the degree by V varies between MINDGV and MAXDGV are removed. */
5771 /* TOUS,AB_SPECIFI:: CARREAU&,CALCUL,&ERREUR */
5773 /* INPUT ARGUMENTS : */
5774 /* ------------------ */
5775 /* NDJACU: Dimension by U of table PATJAC. */
5776 /* NDJACV: Dimension by V of table PATJAC. */
5777 /* NDIMEN: Dimension of the space. */
5778 /* MINDGU: Lower limit of index by U of coeff. of PATJAC to be taken into account. */
5779 /* MAXDGU: Upper limit of index by U of coeff. of PATJAC to be taken into account. */
5780 /* MINDGV: Lower limit of index by V of coeff. of PATJAC to be taken into account. */
5781 /* MAXDGV: Upper limit of index by V of coeff. of PATJAC to be taken into account. */
5782 /* IORDRU: Order of continuity by U provided by square PATJAC (from -1 to 2) */
5783 /* IORDRV: Order of continuity by U provided by square PATJAC (from -1 to 2) */
5784 /* XMAXJU: Maximum value of Jacobi polynoms of order IORDRU, */
5785 /* from degree 0 to MAXDGU - 2*(IORDU+1) */
5786 /* XMAXJV: Maximum value of Jacobi polynoms of order IORDRV, */
5787 /* from degree 0 to MAXDGV - 2*(IORDV+1) */
5788 /* PATJAC: Table of coeff. of square of approximation with */
5789 /* constraints of order IORDRU by U and IORDRV by V. */
5790 /* VECERR: Auxiliary vector. */
5791 /* ERREUR: MAX Error commited during removal of ALREADY CALCULATED coeff of PATJAC */
5793 /* OUTPUT ARGUMENTS : */
5794 /* ------------------- */
5795 /* ERREUR: MAX Error commited during removal of coeff of PATJAC */
5796 /* of indices from MINDGU to MAXDGU by U and from MINDGV to MAXDGV by V */
5797 /* THEN the already calculated error. */
5799 /* COMMONS USED : */
5800 /* ---------------- */
5802 /* REFERENCES CALLED : */
5803 /* --------------------- */
5805 /* DESCRIPTION/NOTES/LIMITATIONS : */
5806 /* ----------------------------------- */
5807 /* Table PATJAC is the place of storage of coeff. Cij of the square of */
5808 /* approximation of F(U,V). The indices i and j show the degree */
5809 /* by U and by V of base polynoms. These polynoms have the form: */
5811 /* ((1 - U*U)**(IORDRU+1)).J(i-2*(IORDRU+1)(U), where */
5813 /* polynom J(i-2*(IORDU+1)(U) is the Jacobi polynom of order */
5814 /* IORDRU+1 (the same by V by replacing U u V in the expression above). */
5816 /* The contribution to the error of term Cij when it is */
5817 /* removed from PATJAC is increased by: */
5819 /* DABS(Cij)*XMAXJU(i-2*(IORDRU+1))*XMAXJV(J-2*(IORDRV+1)) where */
5821 /* XMAXJU(i-2*(IORDRU+1) = ((1 - U*U)**(IORDRU+1)).J(i-2*(IORDRU+1)(U),
5823 /* XMAXJV(i-2*(IORDRV+1) = ((1 - V*V)**(IORDRV+1)).J(j-2*(IORDRV+1)(V).
5827 /* ***********************************************************************
5829 /* Name of the routine */
5832 /* ----------------------------- Initialisations ------------------------
5835 /* Parameter adjustments */
5837 patjac_dim1 = *ndjacu + 1;
5838 patjac_dim2 = *ndjacv + 1;
5839 patjac_offset = patjac_dim1 * patjac_dim2;
5840 patjac -= patjac_offset;
5843 ldbg = AdvApp2Var_SysBase::mnfndeb_() >= 3;
5845 AdvApp2Var_SysBase::mgenmsg_("MMA2ER1", 7L);
5848 minu = (*iordru + 1) << 1;
5849 minv = (*iordrv + 1) << 1;
5851 /* ------------------- Calculate the increment of the max error --------------- */
5852 /* ----- during the removal of the coeffs of indices from MINDGU to MAXDGU ---- */
5853 /* ---------------- by U and indices from MINDGV to MAXDGV by V --------------- */
5856 for (nd = 1; nd <= i__1; ++nd) {
5859 for (jj = *mindgv; jj <= i__2; ++jj) {
5862 for (ii = *mindgu; ii <= i__3; ++ii) {
5863 bid0 += (d__1 = patjac[ii + (jj + nd * patjac_dim2) *
5864 patjac_dim1], advapp_abs(d__1)) * xmaxju[ii - minu];
5867 bid1 = bid0 * xmaxjv[jj - minv] + bid1;
5875 /* ----------------------- Calculate the max error ----------------------*/
5877 bid1 = AdvApp2Var_MathBase::mzsnorm_(ndimen, &vecerr[1]);
5881 *erreur = AdvApp2Var_MathBase::mzsnorm_(&nd, vaux);
5883 /* ------------------------- The end ------------------------------------
5887 AdvApp2Var_SysBase::mgsomsg_("MMA2ER1", 7L);
5892 //=======================================================================
5893 //function : mma2er2_
5895 //=======================================================================
5896 int mma2er2_(integer *ndjacu,
5908 doublereal *epmscut,
5915 /* System generated locals */
5916 integer patjac_dim1, patjac_dim2, patjac_offset, i__1, i__2;
5919 /* Local variables */
5922 integer i2rdu, i2rdv;
5923 doublereal errnu, errnv;
5924 integer ii, nd, jj, nu, nv;
5925 doublereal bid0, bid1;
5927 /* **********************************************************************
5932 /* Remove coefficients of PATJAC to obtain the minimum degree */
5933 /* by U and V checking the imposed tolerance. */
5937 /* TOUS,AB_SPECIFI:: CARREAU&,CALCUL,&ERREUR */
5939 /* INPUT ARGUMENTS : */
5940 /* ------------------ */
5941 /* NDJACU: Degree by U of table PATJAC. */
5942 /* NDJACV: Degree by V of table PATJAC. */
5943 /* NDIMEN: Dimension of the space. */
5944 /* MINDGU: Limit of index by U of coeff. of PATJAC to be PRESERVED (should be >=0). */
5945 /* MAXDGU: Upper limit of index by U of coeff. of PATJAC to be taken into account. */
5946 /* MINDGV: Limit of index by V of coeff. of PATJAC to be PRESERVED (should be >=0). */
5947 /* MAXDGV: Upper limit of index by V of coeff. of PATJAC to be taken into account. */
5948 /* IORDRU: Order of continuity by U provided by square PATJAC (from -1 to 2) */
5949 /* IORDRV: Order of continuity by U provided by square PATJAC (from -1 to 2) */
5950 /* XMAXJU: Maximum value of Jacobi polynoms of order IORDRU, */
5951 /* from degree 0 to MAXDGU - 2*(IORDU+1) */
5952 /* XMAXJV: Maximum value of Jacobi polynoms of order IORDRV, */
5953 /* from degree 0 to MAXDGV - 2*(IORDV+1) */
5954 /* PATJAC: Table of coeff. of square of approximation with */
5955 /* constraints of order IORDRU by U and IORDRV by V. */
5956 /* EPMSCUT: Tolerance of approximation. */
5957 /* VECERR: Auxiliary vector. */
5958 /* ERREUR: MAX Error commited ALREADY CALCULATED */
5960 /* OUTPUT ARGUMENTS : */
5961 /* ------------------- */
5962 /* ERREUR: MAX Error commited by preserving only coeff of PATJAC */
5963 /* of indices from 0 to NEWDGU by U and from 0 to NEWDGV by V */
5964 /* PLUS the already calculated error. */
5965 /* NEWDGU: Min. Degree by U such as the square of approximation */
5966 /* could check the tolerance. There is always NEWDGU >= MINDGU >= 0. */
5967 /* NEWDGV: Min. Degree by V such as the square of approximation */
5968 /* could check the tolerance. There is always NEWDGV >= MINDGV >= 0. */
5971 /* COMMONS USED : */
5972 /* ---------------- */
5974 /* REFERENCES CALLED : */
5975 /* --------------------- */
5977 /* DESCRIPTION/NOTES/LIMITATIONS : */
5978 /* ----------------------------------- */
5979 /* Table PATJAC is the place of storage of coeff. Cij of the square of */
5980 /* approximation of F(U,V). The indices i and j show the degree */
5981 /* by U and by V of base polynoms. These polynoms have the form: */
5983 /* ((1 - U*U)**(IORDRU+1)).J(i-2*(IORDRU+1)(U), where */
5985 /* polynom J(i-2*(IORDU+1)(U) is the Jacobi polynom of order */
5986 /* IORDRU+1 (the same by V by replacing U u V in the expression above). */
5988 /* The contribution to the error of term Cij when it is */
5989 /* removed from PATJAC is increased by: */
5991 /* DABS(Cij)*XMAXJU(i-2*(IORDRU+1))*XMAXJV(J-2*(IORDRV+1)) where */
5993 /* XMAXJU(i-2*(IORDRU+1) = ((1 - U*U)**(IORDRU+1)).J(i-2*(IORDRU+1)(U),
5995 /* XMAXJV(i-2*(IORDRV+1) = ((1 - V*V)**(IORDRV+1)).J(j-2*(IORDRV+1)(V).
5999 /* **********************************************************************
6001 /* Name of the routine */
6004 /* ----------------------------- Initialisations ------------------------
6007 /* Parameter adjustments */
6009 patjac_dim1 = *ndjacu + 1;
6010 patjac_dim2 = *ndjacv + 1;
6011 patjac_offset = patjac_dim1 * patjac_dim2;
6012 patjac -= patjac_offset;
6015 ldbg = AdvApp2Var_SysBase::mnfndeb_() >= 3;
6017 AdvApp2Var_SysBase::mgenmsg_("MMA2ER2", 7L);
6020 i2rdu = (*iordru + 1) << 1;
6021 i2rdv = (*iordrv + 1) << 1;
6025 /* **********************************************************************
6027 /* -------------------- Cutting of oefficients ------------------------
6029 /* **********************************************************************
6034 /* ------------------- Calculate the increment of max error --------------- */
6035 /* ----- during the removal of coeff. of indices from MINDGU to MAXDGU ------ */
6036 /* ---------------- by U, the degree by V is fixed to NV -----------------
6041 bid0 = xmaxjv[nv - i2rdv];
6043 for (nd = 1; nd <= i__1; ++nd) {
6046 for (ii = i2rdu; ii <= i__2; ++ii) {
6047 bid1 += (d__1 = patjac[ii + (nv + nd * patjac_dim2) *
6048 patjac_dim1], advapp_abs(d__1)) * xmaxju[ii - i2rdu] * bid0;
6055 vecerr[1] = *epmscut * 2;
6057 errnv = AdvApp2Var_MathBase::mzsnorm_(ndimen, &vecerr[1]);
6059 /* ------------------- Calculate the increment of max error --------------- */
6060 /* ----- during the removal of coeff. of indices from MINDGV to MAXDGV ------ */
6061 /* ---------------- by V, the degree by U is fixed to NU -----------------
6066 bid0 = xmaxju[nu - i2rdu];
6068 for (nd = 1; nd <= i__1; ++nd) {
6071 for (jj = i2rdv; jj <= i__2; ++jj) {
6072 bid1 += (d__1 = patjac[nu + (jj + nd * patjac_dim2) *
6073 patjac_dim1], advapp_abs(d__1)) * xmaxjv[jj - i2rdv] * bid0;
6080 vecerr[1] = *epmscut * 2;
6082 errnu = AdvApp2Var_MathBase::mzsnorm_(ndimen, &vecerr[1]);
6084 /* ----------------------- Calculate the max error ----------------------
6090 errnu = AdvApp2Var_MathBase::mzsnorm_(&nd, vaux);
6092 errnv = AdvApp2Var_MathBase::mzsnorm_(&nd, vaux);
6094 if (errnu > errnv) {
6095 if (errnv < *epmscut) {
6102 if (errnu < *epmscut) {
6112 /* -------------------------- Return the degrees -------------------
6116 *newdgu = advapp_max(nu,1);
6117 *newdgv = advapp_max(nv,1);
6119 /* ----------------------------------- The end --------------------------
6123 AdvApp2Var_SysBase::mgsomsg_("MMA2ER2", 7L);
6128 //=======================================================================
6129 //function : mma2fnc_
6131 //=======================================================================
6132 int AdvApp2Var_ApproxF2var::mma2fnc_(integer *ndimen,
6136 const AdvApp2Var_EvaluatorFunc2Var& foncnp,
6162 /* System generated locals */
6163 integer courbe_dim1, courbe_dim2, courbe_offset, somtab_dim1, somtab_dim2,
6164 somtab_offset, diftab_dim1, diftab_dim2, diftab_offset,
6165 contr1_dim1, contr1_dim2, contr1_offset, contr2_dim1, contr2_dim2,
6166 contr2_offset, errmax_dim1, errmax_offset, errmoy_dim1,
6167 errmoy_offset, i__1;
6170 /* Local variables */
6173 integer ideb1, ibid1, ibid2, ncfja, ndgre, ilong,
6175 doublereal* wrkar = 0;
6178 doublereal uvpav[4] /* was [2][2] */;
6182 doublereal uv11[4] /* was [2][2] */;
6185 integer isz1, isz2, isz3, isz4, isz5;
6186 intptr_t ipt1, ipt2, ipt3, ipt4, ipt5,iptt, jptt;
6188 /* **********************************************************************
6193 /* Approximation of a limit of non polynomial function F(u,v) */
6194 /* (in the space of dimension NDIMEN) by SEVERAL */
6195 /* polynomial curves, by the method of least squares. The parameter of the function is preserved. */
6199 /* TOUS, AB_SPECIFI :: FONCTION&,EXTREMITE&, APPROXIMATION, &COURBE. */
6201 /* INPUT ARGUMENTS : */
6202 /* ----------------- */
6203 /* NDIMEN: Total Dimension of the space (sum of dimensions */
6204 /* of sub-spaces) */
6205 /* NBSESP: Number of "independent" sub-spaces. */
6206 /* NDIMSE: Table of dimensions of sub-spaces. */
6207 /* UVFONC: Limits of the interval (a,b)x(c,d) of definition of the */
6208 /* function to be approached by U (UVFONC(*,1) contains (a,b)) */
6209 /* and by V (UVFONC(*,2) contains (c,d)). */
6210 /* FONCNP: External function of position on the non polynomial function to be approached. */
6211 /* TCONST: Value of isoparameter of F(u,v) to be discretized. */
6212 /* ISOFAV: Type of chosen iso, = 1, shose that discretization is with u */
6213 /* fixed; = 2, shows that v is fixed. */
6214 /* NBROOT: Nb of points of discretisation of the iso, extremities not included. */
6215 /* ROOTLG: Table of roots of the polynom of Legendre defined on */
6216 /* (-1,1), of degree NBROOT. */
6217 /* IORDRE: Order of constraint at the extremities of the limit */
6218 /* -1 = no constraints, */
6219 /* 0 = constraints of passage to limits (i.e. C0), */
6220 /* 1 = C0 + constraints of 1st derivatives (i.e. C1), */
6221 /* 2 = C1 + constraints of 2nd derivatives (i.e. C2). */
6222 /* IDERIV: Order of derivative of the limit. */
6223 /* NDGJAC: Degree of serial development to be used for calculation in */
6224 /* the Jacobi base. */
6225 /* NBCRMX: Max Nb of curves to be created. */
6226 /* NCFLIM: Max Nb of coeff of the polynomial curve */
6227 /* of approximation (should be above or equal to */
6228 /* 2*IORDRE+2 and below or equal to 50). */
6229 /* EPSAPR: Table of required errors of approximation */
6230 /* sub-space by sub-space. */
6232 /* OUTPUT ARGUMENTS : */
6233 /* ------------------- */
6234 /* NCOEFF: Number of significative coeff of calculated curves. */
6235 /* COURBE: Table of coeff. of calculated polynomial curves. */
6236 /* Should be dimensioned in (NCFLIM,NDIMEN,NBCRMX). */
6237 /* These curves are ALWAYS parametrized in (-1,1). */
6238 /* NBCRBE: Nb of calculated curves. */
6239 /* SOMTAB: For F defined on (-1,1) (otherwise rescale the */
6240 /* parameters), this is the table of sums F(u,vj) + F(u,-vj)
6242 /* if ISOFAV = 1 (and IDERIV=0, otherwise the derivatives */
6243 /* by u of order IDERIV are taken) or sumes F(ui,v) + F(-ui,v) if */
6244 /* ISOFAV = 2 (and IDERIV=0, otherwise the derivatives by */
6245 /* v of order IDERIV are taken). */
6246 /* DIFTAB: For F defined on (-1,1) (otherwise rescale the */
6247 /* parameters), this is the table of sums F(u,vj) - F(u,-vj)
6249 /* if ISOFAV = 1 (and IDERIV=0, otherwise the derivatives */
6250 /* by u of order IDERIV are taken) or sumes F(ui,v) + F(-ui,v) if */
6251 /* ISOFAV = 2 (and IDERIV=0, otherwise the derivatives by */
6252 /* v of order IDERIV are taken). */
6253 /* CONTR1: Contains the coordinates of the left extremity of the iso */
6254 /* and of its derivatives till order IORDRE */
6255 /* CONTR2: Contains the coordinates of the right extremity of the iso */
6256 /* and of its derivatives till order IORDRE */
6257 /* TABDEC: Table of NBCRBE+1 parameters of cut of UVFONC(1:2,1)
6259 /* if ISOFAV=2, or of UVFONC(1:2,2) if ISOFAV=1. */
6260 /* ERRMAX: Table of MAX errors (sub-space by sub-space) */
6261 /* committed in the approximation of FONCNP by NBCRBE curves. */
6262 /* ERRMOY: Table of AVERAGE errors (sub-space by sub-space) */
6263 /* committed in the approximation of FONCNP by NBCRBE curves.
6264 /* IERCOD: Error code: */
6265 /* -1 = ERRMAX > EPSAPR for at least one sub-space. */
6266 /* (the resulting curves of at least mathematic degree NCFLIM-1 */
6267 /* are calculated). */
6268 /* 0 = Everything is ok. */
6269 /* 1 = Pb of incoherence of inputs. */
6270 /* 10 = Pb of calculation of the interpolation of constraints. */
6271 /* 13 = Pb in the dynamic allocation. */
6272 /* 33 = Pb in the data recuperation from block data */
6273 /* of coeff. of integration by GAUSS method. */
6274 /* >100 Pb in the evaluation of FONCNP, the returned error code */
6275 /* is equal to the error code of FONCNP + 100. */
6277 /* COMMONS USED : */
6278 /* ---------------- */
6280 /* REFERENCES CALLED : */
6281 /* ----------------------- */
6283 /* DESCRIPTION/NOTES/LIMITATIONS : */
6284 /* ----------------------------------- */
6285 /* --> The approximation part is done in the space of dimension */
6286 /* NDIMEN (the sum of dimensions of sub-spaces). For example : */
6287 /* If NBSESP=2 and NDIMSE(1)=3, NDIMSE(2)=2, there is smoothing with */
6288 /* NDIMEN=5. The result (in COURBE(NDIMEN,NCOEFF,i) ), will be */
6289 /* composed of the result of smoothing of 3D function in */
6290 /* COURBE(1:3,1:NCOEFF,i) and of smoothing of 2D function in */
6291 /* COURBE(4:5,1:NCOEFF,i). */
6293 /* --> Routine FONCNP should be declared EXTERNAL in the program */
6294 /* calling MMA2FNC. */
6296 /* --> Function FONCNP, declared externally, should be declared */
6297 /* IMPERATIVELY in form : */
6298 /* SUBROUTINE FONCNP(NDIMEN,UINTFN,VINTFN,IIUOUV,TCONST,NBPTAB */
6299 /* ,TTABLE,IDERIU,IDERIV,IERCOD) */
6300 /* where the input arguments are : */
6301 /* - NDIMEN is integer defined as the sum of dimensions of */
6302 /* sub-spaces (i.e. total dimension of the problem). */
6303 /* - UINTFN(2) is a table of 2 reals containing the interval */
6304 /* by u where the function to be approximated is defined */
6305 /* (so it is equal to UIFONC). */
6306 /* - VINTFN(2) is a table of 2 reals containing the interval */
6307 /* by v where the function to be approximated is defined */
6308 /* (so it is equal to VIFONC). */
6309 /* - IIUOUV, shows that the points to be calculated have a constant U */
6310 /* (IIUOUV=1) or a constant V (IIUOUV=2). */
6311 /* - TCONST, real, value of the fixed discretisation parameter. Takes values */
6312 /* in (UINTFN(1),UINTFN(2)) if IIUOUV=1, */
6313 /* or in (VINTFN(1),VINTFN(2)) if IIUOUV=2. */
6314 /* - NBPTAB, the nb of point of discretisation following the free variable */
6315 /* : V if IIUOUV=1 or U if IIUOUV = 2. */
6316 /* - TTABLE, Table of NBPTAB parametres of discretisation. . */
6317 /* - IDERIU, integer, takes values between 0 (position) */
6318 /* and IORDREU (partial derivative of the function by u */
6319 /* of order IORDREU if IORDREU > 0). */
6320 /* - IDERIV, integer, takes values between 0 (position) */
6321 /* and IORDREV (partial derivative of the function by v */
6322 /* of order IORDREV if IORDREV > 0). */
6323 /* and the output arguments are : */
6324 /* - FPNTAB(NDIMEN,NBPTAB) contains, at output, the table of */
6325 /* NBPTAB points calculated in FONCNP. */
6326 /* - IERCOD is, at output the error code of FONCNP. This code */
6327 /* (integer) should be strictly positive if there is a problem. */
6329 /* The input arguments SHOULD NOT BE modified under FONCNP.
6332 /* --> If IERCOD=-1, the required precision can't be reached (ERRMAX */
6333 /* is above EPSAPR on at least one sub-space), but
6335 /* one gives the best possible result for NCFLIM and EPSAPR */
6336 /* chosen by the user. In this case (and for IERCOD=0), there is a solution. */
6339 /* **********************************************************************
6341 /* Name of the routine */
6343 /* Parameter adjustments */
6348 errmoy_dim1 = *nbsesp;
6349 errmoy_offset = errmoy_dim1 + 1;
6350 errmoy -= errmoy_offset;
6351 errmax_dim1 = *nbsesp;
6352 errmax_offset = errmax_dim1 + 1;
6353 errmax -= errmax_offset;
6354 contr2_dim1 = *ndimen;
6355 contr2_dim2 = *iordre + 2;
6356 contr2_offset = contr2_dim1 * (contr2_dim2 + 1) + 1;
6357 contr2 -= contr2_offset;
6358 contr1_dim1 = *ndimen;
6359 contr1_dim2 = *iordre + 2;
6360 contr1_offset = contr1_dim1 * (contr1_dim2 + 1) + 1;
6361 contr1 -= contr1_offset;
6362 diftab_dim1 = *nbroot / 2 + 1;
6363 diftab_dim2 = *ndimen;
6364 diftab_offset = diftab_dim1 * (diftab_dim2 + 1);
6365 diftab -= diftab_offset;
6366 somtab_dim1 = *nbroot / 2 + 1;
6367 somtab_dim2 = *ndimen;
6368 somtab_offset = somtab_dim1 * (somtab_dim2 + 1);
6369 somtab -= somtab_offset;
6371 courbe_dim1 = *ncflim;
6372 courbe_dim2 = *ndimen;
6373 courbe_offset = courbe_dim1 * (courbe_dim2 + 1) + 1;
6374 courbe -= courbe_offset;
6375 AdvApp2Var_SysBase anAdvApp2Var_SysBase;
6378 ibb = AdvApp2Var_SysBase::mnfndeb_();
6380 AdvApp2Var_SysBase::mgenmsg_("MMA2FNC", 7L);
6385 /* ---------------- Set to zero the coefficients of CURVE --------------
6388 ilong = *ndimen * *ncflim * *nbcrmx;
6389 AdvApp2Var_SysBase::mvriraz_(&ilong, &courbe[courbe_offset]);
6391 /* **********************************************************************
6393 /* -------------------------- Checking of entries ------------------
6395 /* **********************************************************************
6398 AdvApp2Var_MathBase::mmveps3_(&eps3);
6399 if ((d__1 = uvfonc[4] - uvfonc[3], advapp_abs(d__1)) < eps3) {
6402 if ((d__1 = uvfonc[6] - uvfonc[5], advapp_abs(d__1)) < eps3) {
6411 /* ********************************************************************** */
6412 /* ------------- Preparation of parameters of discretisation ----------- */
6413 /* **********************************************************************
6416 /* -- Allocation of a table of parameters and points of discretisation -- */
6417 /* --> For the parameters of discretisation. */
6419 /* --> For the points of discretisation in MMA1FDI and MMA1CDI and
6421 /* the auxiliary curve for MMAPCMP */
6422 ibid1 = *ndimen * (*nbroot + 2);
6423 ibid2 = ((*iordre + 1) << 1) * *nbroot;
6424 isz2 = advapp_max(ibid1,ibid2);
6425 ibid1 = (((*ncflim - 1) / 2 + 1) << 1) * *ndimen;
6426 isz2 = advapp_max(ibid1,isz2);
6427 /* --> To return the polynoms of hermit. */
6428 isz3 = ((*iordre + 1) << 2) * (*iordre + 1);
6429 /* --> For the Gauss coeff. of integration. */
6430 isz4 = (*nbroot / 2 + 1) * (*ndgjac + 1 - ((*iordre + 1) << 1));
6431 /* --> For the coeff of the curve in the base of Jacobi */
6432 isz5 = (*ndgjac + 1) * *ndimen;
6434 ndwrk = isz1 + isz2 + isz3 + isz4 + isz5;
6435 anAdvApp2Var_SysBase.mcrrqst_(&c__8, &ndwrk, wrkar, &iofwr, &ier);
6438 /* --> For the parameters of discretisation (NBROOT+2 extremities). */
6440 /* --> For the points of discretisation FPNTAB(NDIMEN,NBROOT+2), */
6441 /* FPNTAB(NBROOT,2*(IORDRE+1)) and for WRKAR of MMAPCMP. */
6443 /* --> For the polynoms of Hermit */
6445 /* --> For the Gauss coeff of integration. */
6447 /* --> For the curve in Jacobi. */
6450 /* ------------------ Initialisation of management of cuts ---------
6454 uvpav[0] = uvfonc[3];
6455 uvpav[1] = uvfonc[4];
6456 tabdec[0] = uvfonc[5];
6457 tabdec[1] = uvfonc[6];
6458 } else if (*isofav == 2) {
6459 tabdec[0] = uvfonc[3];
6460 tabdec[1] = uvfonc[4];
6461 uvpav[2] = uvfonc[5];
6462 uvpav[3] = uvfonc[6];
6470 /* **********************************************************************
6472 /* APPROXIMATION WITH CUTS */
6473 /* **********************************************************************
6477 /* --> When the top is reached, this is the end ! */
6478 if (nupil - *nbcrbe == 0) {
6483 uvpav[2] = tabdec[*nbcrbe];
6484 uvpav[3] = tabdec[*nbcrbe + 1];
6485 } else if (*isofav == 2) {
6486 uvpav[0] = tabdec[*nbcrbe];
6487 uvpav[1] = tabdec[*nbcrbe + 1];
6492 /* -------------------- Normalization of parameters -------------------- */
6494 mma1nop_(nbroot, &rootlg[1], uvpav, isofav, &wrkar[ipt1], &ier);
6499 /* -------------------- Discretisation of FONCNP ------------------------ */
6501 mma1fdi_(ndimen, uvpav, foncnp, isofav, tconst, nbroot, &wrkar[ipt1],
6502 iordre, ideriv, &wrkar[ipt2], &somtab[(ncb1 * somtab_dim2 + 1) *
6503 somtab_dim1], &diftab[(ncb1 * diftab_dim2 + 1) * diftab_dim1], &
6504 contr1[(ncb1 * contr1_dim2 + 1) * contr1_dim1 + 1], &contr2[(ncb1
6505 * contr2_dim2 + 1) * contr2_dim1 + 1], iercod);
6510 /* -----------Cut the discretisation of constraints ------------*/
6513 mma1cdi_(ndimen, nbroot, &rootlg[1], iordre, &contr1[(ncb1 *
6514 contr1_dim2 + 1) * contr1_dim1 + 1], &contr2[(ncb1 *
6515 contr2_dim2 + 1) * contr2_dim1 + 1], &somtab[(ncb1 *
6516 somtab_dim2 + 1) * somtab_dim1], &diftab[(ncb1 * diftab_dim2
6517 + 1) * diftab_dim1], &wrkar[ipt2], &wrkar[ipt3], &ier);
6523 /* **********************************************************************
6525 /* -------------------- Calculate the curve of approximation -------------
6527 /* **********************************************************************
6530 mma1jak_(ndimen, nbroot, iordre, ndgjac, &somtab[(ncb1 * somtab_dim2 + 1)
6531 * somtab_dim1], &diftab[(ncb1 * diftab_dim2 + 1) * diftab_dim1], &
6532 wrkar[ipt4], &wrkar[ipt5], &ier);
6537 /* **********************************************************************
6539 /* ---------------- Add polynom of interpolation -------------------
6541 /* **********************************************************************
6545 mma1cnt_(ndimen, iordre, &contr1[(ncb1 * contr1_dim2 + 1) *
6546 contr1_dim1 + 1], &contr2[(ncb1 * contr2_dim2 + 1) *
6547 contr2_dim1 + 1], &wrkar[ipt3], ndgjac, &wrkar[ipt5]);
6550 /* **********************************************************************
6552 /* --------------- Calculate Max and Average error ----------------------
6554 /* **********************************************************************
6557 mma1fer_(ndimen, nbsesp, &ndimse[1], iordre, ndgjac, &wrkar[ipt5], ncflim,
6558 &epsapr[1], &wrkar[ipt2], &errmax[ncb1 * errmax_dim1 + 1], &
6559 errmoy[ncb1 * errmoy_dim1 + 1], &ncoeff[ncb1], &ier);
6564 if (ier == 0 || (ier == -1 && nupil == *nbcrmx)) {
6566 /* ******************************************************************
6568 /* ----------------------- Compression du resultat ------------------
6570 /* ******************************************************************
6576 ncfja = *ndgjac + 1;
6577 /* -> Compression of result in WRKAR(IPT2) */
6580 AdvApp2Var_MathBase::mmapcmp_(ndimen,
6581 &ncfja, &ncoeff[ncb1], &wrkar[ipt5], &wrkar[ipt2]);
6583 AdvApp2Var_MathBase::mmapcmp_((integer*)ndimen,
6589 ilong = *ndimen * *ncflim;
6590 AdvApp2Var_SysBase::mvriraz_(&ilong, &wrkar[ipt5]);
6591 /* -> Passage to canonic base (-1,1) (result in WRKAR(IPT5)).
6593 ndgre = ncoeff[ncb1] - 1;
6595 for (nd = 1; nd <= i__1; ++nd) {
6596 iptt = ipt2 + ((nd - 1) << 1) * (ndgre / 2 + 1);
6597 jptt = ipt5 + (nd - 1) * ncoeff[ncb1];
6598 AdvApp2Var_MathBase::mmjacan_(iordre, &ndgre, &wrkar[iptt], &wrkar[jptt]);
6602 /* -> Store the calculated curve */
6604 AdvApp2Var_MathBase::mmfmca8_(&ncoeff[ncb1], ndimen, &ibid1, ncflim, ndimen, &ibid1, &
6605 wrkar[ipt5], &courbe[(ncb1 * courbe_dim2 + 1) * courbe_dim1 +
6608 /* -> Before normalization of constraints on (-1,1), recalculate */
6609 /* the true constraints. */
6611 for (ii = 0; ii <= i__1; ++ii) {
6612 mma1noc_(uv11, ndimen, &ii, &contr1[(ii + 1 + ncb1 * contr1_dim2)
6613 * contr1_dim1 + 1], uvpav, isofav, ideriv, &contr1[(ii +
6614 1 + ncb1 * contr1_dim2) * contr1_dim1 + 1]);
6615 mma1noc_(uv11, ndimen, &ii, &contr2[(ii + 1 + ncb1 * contr2_dim2)
6616 * contr2_dim1 + 1], uvpav, isofav, ideriv, &contr2[(ii +
6617 1 + ncb1 * contr2_dim2) * contr2_dim1 + 1]);
6621 ibid1 = (*nbroot / 2 + 1) * *ndimen;
6622 mma1noc_(uv11, &ibid1, &ii, &somtab[(ncb1 * somtab_dim2 + 1) *
6623 somtab_dim1], uvpav, isofav, ideriv, &somtab[(ncb1 *
6624 somtab_dim2 + 1) * somtab_dim1]);
6625 mma1noc_(uv11, &ibid1, &ii, &diftab[(ncb1 * diftab_dim2 + 1) *
6626 diftab_dim1], uvpav, isofav, ideriv, &diftab[(ncb1 *
6627 diftab_dim2 + 1) * diftab_dim1]);
6630 for (nd = 1; nd <= i__1; ++nd) {
6631 mma1noc_(uv11, &ncoeff[ncb1], &ii, &courbe[(nd + ncb1 *
6632 courbe_dim2) * courbe_dim1 + 1], uvpav, isofav, ideriv, &
6633 courbe[(nd + ncb1 * courbe_dim2) * courbe_dim1 + 1]);
6637 /* -> Update the nb of already created curves */
6640 /* -> ...otherwise try to cut the current interval in 2... */
6642 tmil = (tabdec[*nbcrbe + 1] + tabdec[*nbcrbe]) / 2.;
6645 ilong = (nupil - *nbcrbe) << 3;
6646 AdvApp2Var_SysBase::mcrfill_(&ilong, &tabdec[ideb],&tabdec[ideb1]);
6647 tabdec[ideb] = tmil;
6651 /* ---------- Make approximation of the rest -----------
6656 /* --------------------- Return code of error -----------------
6658 /* --> Pb with dynamic allocation */
6662 /* --> Inputs incoherent. */
6667 /* -------------------------- Dynamic desallocation -------------------
6672 anAdvApp2Var_SysBase.mcrdelt_(&c__8, &ndwrk, wrkar, &iofwr, &ier);
6679 /* ------------------------------ The end -------------------------------
6684 AdvApp2Var_SysBase::maermsg_("MMA2FNC", iercod, 7L);
6687 AdvApp2Var_SysBase::mgsomsg_("MMA2FNC", 7L);
6692 //=======================================================================
6693 //function : mma2fx6_
6695 //=======================================================================
6696 int AdvApp2Var_ApproxF2var::mma2fx6_(integer *ncfmxu,
6713 /* System generated locals */
6714 integer epsfro_dim1, epsfro_offset, patcan_dim1, patcan_dim2, patcan_dim3,
6715 patcan_dim4, patcan_offset, errmax_dim1, errmax_dim2,
6716 errmax_offset, ncoefu_dim1, ncoefu_offset, ncoefv_dim1,
6717 ncoefv_offset, i__1, i__2, i__3, i__4, i__5;
6718 doublereal d__1, d__2;
6720 /* Local variables */
6721 integer idim, ncfu, ncfv, id, ii, nd, jj, ku, kv, ns, ibb;
6725 /* **********************************************************************
6730 /* Reduction of degree when the squares are the squares of constraints. */
6734 /* TOUS,AB_SPECIFI::CARREAU&,REDUCTION,&CARREAU */
6736 /* INPUT ARGUMENTS : */
6737 /* ------------------ */
6738 /* NCFMXU: Max Nb of coeff by u of solution P(u,v) (table */
6739 /* PATCAN). This argument serves only to declare the size of this table. */
6740 /* NCFMXV: Max Nb of coeff by v of solution P(u,v) (table */
6741 /* PATCAN). This argument serves only to declare the size of this table. */
6742 /* NDIMEN: Total dimension of the space where the processed function */
6743 /* takes its values.(sum of dimensions of sub-spaces) */
6744 /* NBSESP: Nb of independent sub-spaces where the errors are measured. */
6745 /* NDIMSE: Table of dimensions of NBSESP sub-spaces. */
6746 /* NBUPAT: Nb of square solution by u. */
6747 /* NBVPAT: Nb of square solution by v. */
6748 /* IORDRU: Order of constraint imposed at the extremities of iso-V */
6749 /* = 0, the extremities of iso-V are calculated */
6750 /* = 1, additionally the 1st derivative in the direction of iso-V is calculated */
6751 /* = 2, additionally the 2nd derivative in the direction of iso-V is calculated */
6752 /* IORDRV: Ordre de contrainte impose aux extremites de l'iso-U */
6753 /* = 0, on calcule les extremites de l'iso-U. */
6754 /* = 1, additionally the 1st derivative in the direction of iso-U is calculated */
6755 /* = 2, additionally the 2nd derivative in the direction of iso-U is calculated */
6756 /* EPSAPR: Table of imposed precisions, sub-space by sub-space. */
6757 /* EPSFRO: Table of imposed precisions, sub-space by sub-space on the limits of squares. */
6758 /* PATCAN: Table of coeff. in the canonic base of squares P(u,v) calculated for (u,v) in (-1,1). */
6759 /* ERRMAX: Table of MAX errors (sub-space by sub-space) */
6760 /* committed in the approximation of F(u,v) by P(u,v). */
6761 /* NCOEFU: Table of Nb of significative coeffs. by u of calculated squares. */
6762 /* NCOEFV: Table of Nb of significative coeffs. by v of calculated squares. */
6764 /* OUTPUT ARGUMENTS : */
6765 /* ------------------- */
6766 /* NCOEFU: Table of Nb of significative coeffs. by u of calculated squares. */
6767 /* NCOEFV: Table of Nb of significative coeffs. by v of calculated squares. */
6769 /* COMMONS USED : */
6770 /* ---------------- */
6772 /* REFERENCES CALLED : */
6773 /* --------------------- */
6775 /* DESCRIPTION/NOTES/LIMITATIONS : */
6776 /* ------------------------------- */
6778 /* **********************************************************************
6781 /* Name of the routine */
6784 /* Parameter adjustments */
6785 epsfro_dim1 = *nbsesp;
6786 epsfro_offset = epsfro_dim1 * 5 + 1;
6787 epsfro -= epsfro_offset;
6790 ncoefv_dim1 = *nbupat;
6791 ncoefv_offset = ncoefv_dim1 + 1;
6792 ncoefv -= ncoefv_offset;
6793 ncoefu_dim1 = *nbupat;
6794 ncoefu_offset = ncoefu_dim1 + 1;
6795 ncoefu -= ncoefu_offset;
6796 errmax_dim1 = *nbsesp;
6797 errmax_dim2 = *nbupat;
6798 errmax_offset = errmax_dim1 * (errmax_dim2 + 1) + 1;
6799 errmax -= errmax_offset;
6800 patcan_dim1 = *ncfmxu;
6801 patcan_dim2 = *ncfmxv;
6802 patcan_dim3 = *ndimen;
6803 patcan_dim4 = *nbupat;
6804 patcan_offset = patcan_dim1 * (patcan_dim2 * (patcan_dim3 * (patcan_dim4
6806 patcan -= patcan_offset;
6809 ibb = AdvApp2Var_SysBase::mnfndeb_();
6811 AdvApp2Var_SysBase::mgenmsg_("MMA2FX6", 7L);
6816 for (jj = 1; jj <= i__1; ++jj) {
6818 for (ii = 1; ii <= i__2; ++ii) {
6819 ncfu = ncoefu[ii + jj * ncoefu_dim1];
6820 ncfv = ncoefv[ii + jj * ncoefv_dim1];
6822 /* ********************************************************************** */
6823 /* -------------------- Reduction of degree by U ------------------------- */
6824 /* ********************************************************************** */
6827 if (ncfu <= (*iordru + 1) << 1 && ncfu > 2) {
6831 for (ns = 1; ns <= i__3; ++ns) {
6834 d__1 = tol, d__2 = epsfro[ns + epsfro_dim1 * 9];
6835 tol = advapp_min(d__1,d__2);
6837 d__1 = tol, d__2 = epsfro[ns + epsfro_dim1 * 10];
6838 tol = advapp_min(d__1,d__2);
6840 d__1 = tol, d__2 = epsfro[ns + epsfro_dim1 * 11];
6841 tol = advapp_min(d__1,d__2);
6843 d__1 = tol, d__2 = epsfro[ns + epsfro_dim1 * 12];
6844 tol = advapp_min(d__1,d__2);
6845 if (ii == 1 || ii == *nbupat || jj == 1 || jj == *nbvpat)
6848 d__1 = tol, d__2 = epsfro[ns + epsfro_dim1 * 5];
6849 tol = advapp_min(d__1,d__2);
6851 d__1 = tol, d__2 = epsfro[ns + epsfro_dim1 * 6];
6852 tol = advapp_min(d__1,d__2);
6854 d__1 = tol, d__2 = epsfro[ns + epsfro_dim1 * 7];
6855 tol = advapp_min(d__1,d__2);
6857 d__1 = tol, d__2 = epsfro[ns + (epsfro_dim1 << 3)];
6858 tol = advapp_min(d__1,d__2);
6863 for (nd = 1; nd <= i__4; ++nd) {
6866 for (kv = 1; kv <= i__5; ++kv) {
6867 bid += (d__1 = patcan[ncfu + (kv + (id + (ii + jj
6868 * patcan_dim4) * patcan_dim3) *
6869 patcan_dim2) * patcan_dim1], advapp_abs(d__1));
6875 if (bid > tol * 1e-6 || bid > errmax[ns + (ii + jj *
6876 errmax_dim2) * errmax_dim1]) {
6887 /* ********************************************************************** */
6888 /* -------------------- Reduction of degree by V ------------------------- */
6889 /* ********************************************************************** */
6892 if (ncfv <= (*iordrv + 1) << 1 && ncfv > 2) {
6896 for (ns = 1; ns <= i__3; ++ns) {
6899 d__1 = tol, d__2 = epsfro[ns + epsfro_dim1 * 9];
6900 tol = advapp_min(d__1,d__2);
6902 d__1 = tol, d__2 = epsfro[ns + epsfro_dim1 * 10];
6903 tol = advapp_min(d__1,d__2);
6905 d__1 = tol, d__2 = epsfro[ns + epsfro_dim1 * 11];
6906 tol = advapp_min(d__1,d__2);
6908 d__1 = tol, d__2 = epsfro[ns + epsfro_dim1 * 12];
6909 tol = advapp_min(d__1,d__2);
6910 if (ii == 1 || ii == *nbupat || jj == 1 || jj == *nbvpat)
6913 d__1 = tol, d__2 = epsfro[ns + epsfro_dim1 * 5];
6914 tol = advapp_min(d__1,d__2);
6916 d__1 = tol, d__2 = epsfro[ns + epsfro_dim1 * 6];
6917 tol = advapp_min(d__1,d__2);
6919 d__1 = tol, d__2 = epsfro[ns + epsfro_dim1 * 7];
6920 tol = advapp_min(d__1,d__2);
6922 d__1 = tol, d__2 = epsfro[ns + (epsfro_dim1 << 3)];
6923 tol = advapp_min(d__1,d__2);
6928 for (nd = 1; nd <= i__4; ++nd) {
6931 for (ku = 1; ku <= i__5; ++ku) {
6932 bid += (d__1 = patcan[ku + (ncfv + (id + (ii + jj
6933 * patcan_dim4) * patcan_dim3) *
6934 patcan_dim2) * patcan_dim1], advapp_abs(d__1));
6940 if (bid > tol * 1e-6 || bid > errmax[ns + (ii + jj *
6941 errmax_dim2) * errmax_dim1]) {
6952 /* --- Return the nbs of coeff. and pass to the next square --- */
6955 ncoefu[ii + jj * ncoefu_dim1] = advapp_max(ncfu,2);
6956 ncoefv[ii + jj * ncoefv_dim1] = advapp_max(ncfv,2);
6962 /* ------------------------------ The End -------------------------------
6966 AdvApp2Var_SysBase::mgsomsg_("MMA2FX6", 7L);
6972 //=======================================================================
6973 //function : mma2jmx_
6975 //=======================================================================
6976 int AdvApp2Var_ApproxF2var::mma2jmx_(integer *ndgjac,
6980 /* Initialized data */
6982 static doublereal xmax2[57] = { .9682458365518542212948163499456,
6983 .986013297183269340427888048593603,
6984 1.07810420343739860362585159028115,
6985 1.17325804490920057010925920756025,
6986 1.26476561266905634732910520370741,
6987 1.35169950227289626684434056681946,
6988 1.43424378958284137759129885012494,
6989 1.51281316274895465689402798226634,
6990 1.5878364329591908800533936587012,
6991 1.65970112228228167018443636171226,
6992 1.72874345388622461848433443013543,
6993 1.7952515611463877544077632304216,
6994 1.85947199025328260370244491818047,
6995 1.92161634324190018916351663207101,
6996 1.98186713586472025397859895825157,
6997 2.04038269834980146276967984252188,
6998 2.09730119173852573441223706382076,
6999 2.15274387655763462685970799663412,
7000 2.20681777186342079455059961912859,
7001 2.25961782459354604684402726624239,
7002 2.31122868752403808176824020121524,
7003 2.36172618435386566570998793688131,
7004 2.41117852396114589446497298177554,
7005 2.45964731268663657873849811095449,
7006 2.50718840313973523778244737914028,
7007 2.55385260994795361951813645784034,
7008 2.59968631659221867834697883938297,
7009 2.64473199258285846332860663371298,
7010 2.68902863641518586789566216064557,
7011 2.73261215675199397407027673053895,
7012 2.77551570192374483822124304745691,
7013 2.8177699459714315371037628127545,
7014 2.85940333797200948896046563785957,
7015 2.90044232019793636101516293333324,
7016 2.94091151970640874812265419871976,
7017 2.98083391718088702956696303389061,
7018 3.02023099621926980436221568258656,
7019 3.05912287574998661724731962377847,
7020 3.09752842783622025614245706196447,
7021 3.13546538278134559341444834866301,
7022 3.17295042316122606504398054547289,
7023 3.2099992681699613513775259670214,
7024 3.24662674946606137764916854570219,
7025 3.28284687953866689817670991319787,
7026 3.31867291347259485044591136879087,
7027 3.35411740487202127264475726990106,
7028 3.38919225660177218727305224515862,
7029 3.42390876691942143189170489271753,
7030 3.45827767149820230182596660024454,
7031 3.49230918177808483937957161007792,
7032 3.5260130200285724149540352829756,
7033 3.55939845146044235497103883695448,
7034 3.59247431368364585025958062194665,
7035 3.62524904377393592090180712976368,
7036 3.65773070318071087226169680450936,
7037 3.68992700068237648299565823810245,
7038 3.72184531357268220291630708234186 };
7039 static doublereal xmax4[55] = { 1.1092649593311780079813740546678,
7040 1.05299572648705464724876659688996,
7041 1.0949715351434178709281698645813,
7042 1.15078388379719068145021100764647,
7043 1.2094863084718701596278219811869,
7044 1.26806623151369531323304177532868,
7045 1.32549784426476978866302826176202,
7046 1.38142537365039019558329304432581,
7047 1.43575531950773585146867625840552,
7048 1.48850442653629641402403231015299,
7049 1.53973611681876234549146350844736,
7050 1.58953193485272191557448229046492,
7051 1.63797820416306624705258190017418,
7052 1.68515974143594899185621942934906,
7053 1.73115699602477936547107755854868,
7054 1.77604489805513552087086912113251,
7055 1.81989256661534438347398400420601,
7056 1.86276344480103110090865609776681,
7057 1.90471563564740808542244678597105,
7058 1.94580231994751044968731427898046,
7059 1.98607219357764450634552790950067,
7060 2.02556989246317857340333585562678,
7061 2.06433638992049685189059517340452,
7062 2.10240936014742726236706004607473,
7063 2.13982350649113222745523925190532,
7064 2.17661085564771614285379929798896,
7065 2.21280102016879766322589373557048,
7066 2.2484214321456956597803794333791,
7067 2.28349755104077956674135810027654,
7068 2.31805304852593774867640120860446,
7069 2.35210997297725685169643559615022,
7070 2.38568889602346315560143377261814,
7071 2.41880904328694215730192284109322,
7072 2.45148841120796359750021227795539,
7073 2.48374387161372199992570528025315,
7074 2.5155912654873773953959098501893,
7075 2.54704548720896557684101746505398,
7076 2.57812056037881628390134077704127,
7077 2.60882970619319538196517982945269,
7078 2.63918540521920497868347679257107,
7079 2.66919945330942891495458446613851,
7080 2.69888301230439621709803756505788,
7081 2.72824665609081486737132853370048,
7082 2.75730041251405791603760003778285,
7083 2.78605380158311346185098508516203,
7084 2.81451587035387403267676338931454,
7085 2.84269522483114290814009184272637,
7086 2.87060005919012917988363332454033,
7087 2.89823818258367657739520912946934,
7088 2.92561704377132528239806135133273,
7089 2.95274375377994262301217318010209,
7090 2.97962510678256471794289060402033,
7091 3.00626759936182712291041810228171,
7092 3.03267744830655121818899164295959,
7093 3.05886060707437081434964933864149 };
7094 static doublereal xmax6[53] = { 1.21091229812484768570102219548814,
7095 1.11626917091567929907256116528817,
7096 1.1327140810290884106278510474203,
7097 1.1679452722668028753522098022171,
7098 1.20910611986279066645602153641334,
7099 1.25228283758701572089625983127043,
7100 1.29591971597287895911380446311508,
7101 1.3393138157481884258308028584917,
7102 1.3821288728999671920677617491385,
7103 1.42420414683357356104823573391816,
7104 1.46546895108549501306970087318319,
7105 1.50590085198398789708599726315869,
7106 1.54550385142820987194251585145013,
7107 1.58429644271680300005206185490937,
7108 1.62230484071440103826322971668038,
7109 1.65955905239130512405565733793667,
7110 1.69609056468292429853775667485212,
7111 1.73193098017228915881592458573809,
7112 1.7671112206990325429863426635397,
7113 1.80166107681586964987277458875667,
7114 1.83560897003644959204940535551721,
7115 1.86898184653271388435058371983316,
7116 1.90180515174518670797686768515502,
7117 1.93410285411785808749237200054739,
7118 1.96589749778987993293150856865539,
7119 1.99721027139062501070081653790635,
7120 2.02806108474738744005306947877164,
7121 2.05846864831762572089033752595401,
7122 2.08845055210580131460156962214748,
7123 2.11802334209486194329576724042253,
7124 2.14720259305166593214642386780469,
7125 2.17600297710595096918495785742803,
7126 2.20443832785205516555772788192013,
7127 2.2325216999457379530416998244706,
7128 2.2602654243075083168599953074345,
7129 2.28768115912702794202525264301585,
7130 2.3147799369092684021274946755348,
7131 2.34157220782483457076721300512406,
7132 2.36806787963276257263034969490066,
7133 2.39427635443992520016789041085844,
7134 2.42020656255081863955040620243062,
7135 2.44586699364757383088888037359254,
7136 2.47126572552427660024678584642791,
7137 2.49641045058324178349347438430311,
7138 2.52130850028451113942299097584818,
7139 2.54596686772399937214920135190177,
7140 2.5703922285006754089328998222275,
7141 2.59459096001908861492582631591134,
7142 2.61856915936049852435394597597773,
7143 2.64233265984385295286445444361827,
7144 2.66588704638685848486056711408168,
7145 2.68923766976735295746679957665724,
7146 2.71238965987606292679677228666411 };
7148 /* System generated locals */
7151 /* Local variables */
7157 /* **********************************************************************
7162 /* Calculate the max of Jacobo polynoms multiplied by the weight on */
7163 /* (-1,1) for order 0,4,6 or Legendre. */
7167 /* LEGENDRE,APPROXIMATION,ERREUR. */
7169 /* INPUT ARGUMENTS : */
7170 /* ------------------ */
7171 /* NDGJAC: Nb of Jacobi coeff. of approximation. */
7172 /* IORDRE: Order of continuity (from -1 to 2) */
7174 /* OUTPUT ARGUMENTS : */
7175 /* ------------------- */
7176 /* XJACMX: Table of maximums of Jacobi polynoms. */
7178 /* COMMONS USED : */
7179 /* ---------------- */
7181 /* REFERENCES CALLED : */
7182 /* --------------------- */
7184 /* DESCRIPTION/NOTES/LIMITATIONS : */
7185 /* ----------------------------------- */
7188 /* ***********************************************************************
7190 /* Name of the routine */
7191 /* ----------------------------- Initialisations ------------------------
7194 ldbg = AdvApp2Var_SysBase::mnfndeb_() >= 3;
7196 AdvApp2Var_SysBase::mgenmsg_("MMA2JMX", 7L);
7199 numax = *ndgjac - ((*iordre + 1) << 1);
7200 if (*iordre == -1) {
7202 for (ii = 0; ii <= i__1; ++ii) {
7203 bid = (ii * 2. + 1.) / 2.;
7204 xjacmx[ii] = sqrt(bid);
7207 } else if (*iordre == 0) {
7209 for (ii = 0; ii <= i__1; ++ii) {
7210 xjacmx[ii] = xmax2[ii];
7213 } else if (*iordre == 1) {
7215 for (ii = 0; ii <= i__1; ++ii) {
7216 xjacmx[ii] = xmax4[ii];
7219 } else if (*iordre == 2) {
7221 for (ii = 0; ii <= i__1; ++ii) {
7222 xjacmx[ii] = xmax6[ii];
7227 /* ------------------------- The end ------------------------------------
7231 AdvApp2Var_SysBase::mgsomsg_("MMA2JMX", 7L);
7236 //=======================================================================
7237 //function : mma2moy_
7239 //=======================================================================
7240 int mma2moy_(integer *ndgumx,
7252 /* System generated locals */
7253 integer patjac_dim1, patjac_dim2, patjac_offset, i__1, i__2, i__3;
7255 /* Local variables */
7257 integer minu, minv, idebu, idebv, ii, nd, jj;
7258 doublereal bid0, bid1;
7261 /* **********************************************************************
7266 /* Calculate the average approximation error made when only */
7267 /* the coefficients of PATJAC of degree between */
7268 /* 2*(IORDRU+1) and MINDGU by U and 2*(IORDRV+1) and MINDGV by V are preserved. */
7272 /* LEGENDRE,APPROXIMATION, AVERAGE ERROR */
7274 /* INPUT ARGUMENTS : */
7275 /* ------------------ */
7276 /* NDGUMX: Dimension by U of table PATJAC. */
7277 /* NDGVMX: Dimension by V of table PATJAC. */
7278 /* NDIMEN: Dimension of the space. */
7279 /* MINDGU: Lower limit of the index by U of PATJAC coeff to be taken into account. */
7280 /* MAXDGU: Upper limit of the index by U of PATJAC coeff to be taken into account. */
7281 /* MINDGV: Lower limit of the index by V of PATJAC coeff to be taken into account. */
7282 /* MAXDGV: Upper limit of the index by V of PATJAC coeff to be taken into account. */
7283 /* IORDRU: Order of continuity by U provided by square PATJAC (from -1 to 2) */
7284 /* IORDRV: Order of continuity by V provided by square PATJAC (from -1 to 2) */
7285 /* PATJAC: Table of coeff. of the approximation square with */
7286 /* constraints of order IORDRU by U and IORDRV by V. */
7288 /* OUTPUT ARGUMENTS : */
7289 /* ------------------- */
7290 /* ERRMOY: Average error commited by preserving only the coeff of */
7291 /* PATJAC 2*(IORDRU+1) in MINDGU by U and 2*(IORDRV+1) in MINDGV by V. */
7293 /* COMMONS USED : */
7294 /* ---------------- */
7296 /* REFERENCES CALLED : */
7297 /* --------------------- */
7299 /* DESCRIPTION/NOTES/LIMITATIONS : */
7300 /* ----------------------------------- */
7301 /* Table PATJAC stores the coeff. Cij of */
7302 /* approximation square F(U,V). Indexes i and j show the degree by */
7303 /* U and by V of the base polynoms. These base polynoms are in the form: */
7305 /* ((1 - U*U)**(IORDRU+1)).J(i-2*(IORDRU+1)(U), where */
7307 /* polynom J(i-2*(IORDU+1)(U) is the Jacobi polynom of order */
7308 /* IORDRU+1 (the same by V by replacing U by V in the above expression). */
7310 /* The contribution to the average error of term Cij when */
7311 /* it is removed from PATJAC is Cij*Cij. */
7314 /* ***********************************************************************
7316 /* Name of the routine */
7319 /* ----------------------------- Initialisations ------------------------
7322 /* Parameter adjustments */
7323 patjac_dim1 = *ndgumx + 1;
7324 patjac_dim2 = *ndgvmx + 1;
7325 patjac_offset = patjac_dim1 * patjac_dim2;
7326 patjac -= patjac_offset;
7329 ldbg = AdvApp2Var_SysBase::mnfndeb_() >= 3;
7331 AdvApp2Var_SysBase::mgenmsg_("MMA2MOY", 7L);
7334 idebu = (*iordru + 1) << 1;
7335 idebv = (*iordrv + 1) << 1;
7336 minu = advapp_max(idebu,*mindgu);
7337 minv = advapp_max(idebv,*mindgv);
7341 /* ------------------ Calculation of the upper bound of the average error ------------ */
7342 /* -------------------- when the coeff. of indexes from MINDGU to MAXDGU ------ */
7343 /* ---------------- by U and of indexes from MINDGV to MAXDGV by V are removed -------------- */
7346 for (nd = 1; nd <= i__1; ++nd) {
7348 for (jj = minv; jj <= i__2; ++jj) {
7350 for (ii = idebu; ii <= i__3; ++ii) {
7351 bid1 = patjac[ii + (jj + nd * patjac_dim2) * patjac_dim1];
7352 bid0 += bid1 * bid1;
7361 for (nd = 1; nd <= i__1; ++nd) {
7363 for (jj = idebv; jj <= i__2; ++jj) {
7365 for (ii = minu; ii <= i__3; ++ii) {
7366 bid1 = patjac[ii + (jj + nd * patjac_dim2) * patjac_dim1];
7367 bid0 += bid1 * bid1;
7375 /* ----------------------- Calculation of the average error -------------
7379 *errmoy = sqrt(bid0);
7381 /* ------------------------- The end ------------------------------------
7385 AdvApp2Var_SysBase::mgsomsg_("MMA2MOY", 7L);
7390 //=======================================================================
7391 //function : mma2roo_
7393 //=======================================================================
7394 int AdvApp2Var_ApproxF2var::mma2roo_(integer *nbpntu,
7399 /* System generated locals */
7402 /* Local variables */
7405 /* **********************************************************************
7410 /* Return roots of Legendre for discretisations. */
7414 /* TOUS, AB_SPECIFI::CONTRAINTE&, DISCRETISATION, &POINT */
7416 /* INPUT ARGUMENTS : */
7417 /* ------------------ */
7418 /* NBPNTU: Nb of INTERNAL parameters of discretization BY U. */
7419 /* This is also the nb of root of the Legendre polynom where the discretization is done. */
7420 /* NBPNTV: Nb of INTERNAL parameters of discretization BY V. */
7421 /* This is also the nb of root of the Legendre polynom where the discretization is done. */
7423 /* OUTPUT ARGUMENTS : */
7424 /* ------------------- */
7425 /* UROOTL: Table of parameters of discretisation ON (-1,1) BY U.
7427 /* VROOTL: Table of parameters of discretisation ON (-1,1) BY V.
7430 /* COMMONS USED : */
7431 /* ---------------- */
7433 /* REFERENCES CALLED : */
7434 /* --------------------- */
7436 /* DESCRIPTION/NOTES/LIMITATIONS : */
7437 /* ----------------------------------- */
7440 /* **********************************************************************
7443 /* Name of the routine */
7446 /* Parameter adjustments */
7451 ibb = AdvApp2Var_SysBase::mnfndeb_();
7453 AdvApp2Var_SysBase::mgenmsg_("MMA2ROO", 7L);
7456 /* ---------------- Return the POSITIVE roots on U ------------------
7459 AdvApp2Var_MathBase::mmrtptt_(nbpntu, &urootl[(*nbpntu + 1) / 2 + 1]);
7461 for (ii = 1; ii <= i__1; ++ii) {
7462 urootl[ii] = -urootl[*nbpntu - ii + 1];
7465 if (*nbpntu % 2 == 1) {
7466 urootl[*nbpntu / 2 + 1] = 0.;
7469 /* ---------------- Return the POSITIVE roots on V ------------------
7472 AdvApp2Var_MathBase::mmrtptt_(nbpntv, &vrootl[(*nbpntv + 1) / 2 + 1]);
7474 for (ii = 1; ii <= i__1; ++ii) {
7475 vrootl[ii] = -vrootl[*nbpntv - ii + 1];
7478 if (*nbpntv % 2 == 1) {
7479 vrootl[*nbpntv / 2 + 1] = 0.;
7482 /* ------------------------------ The End -------------------------------
7486 AdvApp2Var_SysBase::mgsomsg_("MMA2ROO", 7L);
7490 //=======================================================================
7491 //function : mmmapcoe_
7493 //=======================================================================
7494 int mmmapcoe_(integer *ndim,
7504 /* System generated locals */
7505 integer somtab_dim1, somtab_offset, diftab_dim1, diftab_offset,
7506 crvjac_dim1, crvjac_offset, gsstab_dim1, i__1, i__2, i__3;
7508 /* Local variables */
7509 integer igss, ikdeb;
7511 integer nd, ik, ir, nbroot, ibb;
7513 /* **********************************************************************
7518 /* Calculate the coefficients of polinomial approximation curve */
7519 /* of degree NDGJAC by the method of smallest squares starting from */
7520 /* the discretization of function on the roots of Legendre polynom */
7521 /* of degree NBPNTS. */
7525 /* FONCTION,APPROXIMATION,COEFFICIENT,POLYNOME */
7527 /* INPUT ARGUMENTS : */
7528 /* ------------------ */
7529 /* NDIM : Dimension of the space. */
7530 /* NDGJAC : Max Degree of the polynom of approximation. */
7531 /* The representation in the orthogonal base starts from degree */
7532 /* 0 to degree NDGJAC-2*(JORDRE+1). The polynomial base */
7533 /* is the base of Jacobi of order -1 (Legendre), 0, 1 and 2 */
7534 /* IORDRE : Order of the base of Jacobi (-1,0,1 or 2). Corresponds */
7535 /* to step of constraints, C0,C1 or C2. */
7536 /* NBPNTS : Degree of the polynom of Legendre on the roots which of */
7537 /* are calculated the coefficients of integration by */
7538 /* Gauss method. It is required to set NBPNTS=30,40,50 or 61 */
7539 /* and NDGJAC < NBPNTS. */
7540 /* SOMTAB : Table of F(ti)+F(-ti) with ti in ROOTAB. */
7541 /* DIFTAB : Table of F(ti)-F(-ti) with ti in ROOTAB. */
7542 /* GSSTAB(i,k) : Table of coefficients of integration by the Gauss method : */
7543 /* i varies from 0 to NBPNTS and */
7544 /* k varies from 0 to NDGJAC-2*(JORDRE+1). */
7546 /* OUTPUT ARGUMENTSE : */
7547 /* ------------------- */
7548 /* CRVJAC : Curve of approximation of FONCNP with eventually */
7549 /* taking into account of constraints at the extremities. */
7550 /* This curve is of degree NDGJAC. */
7552 /* COMMONS USED : */
7553 /* ---------------- */
7555 /* REFERENCES CALLED : */
7556 /* --------------------- */
7558 /* DESCRIPTION/NOTES/LIMITATIONS : */
7559 /* ------------------------------- */
7561 /* **********************************************************************
7564 /* Name of the routine */
7566 /* Parameter adjustments */
7567 crvjac_dim1 = *ndgjac + 1;
7568 crvjac_offset = crvjac_dim1;
7569 crvjac -= crvjac_offset;
7570 gsstab_dim1 = *nbpnts / 2 + 1;
7571 diftab_dim1 = *nbpnts / 2 + 1;
7572 diftab_offset = diftab_dim1;
7573 diftab -= diftab_offset;
7574 somtab_dim1 = *nbpnts / 2 + 1;
7575 somtab_offset = somtab_dim1;
7576 somtab -= somtab_offset;
7579 ibb = AdvApp2Var_SysBase::mnfndeb_();
7581 AdvApp2Var_SysBase::mgenmsg_("MMMAPCO", 7L);
7583 ikdeb = (*iordre + 1) << 1;
7584 nbroot = *nbpnts / 2;
7587 for (nd = 1; nd <= i__1; ++nd) {
7589 /* ----------------- Calculate the coefficients of even degree ----------
7593 for (ik = ikdeb; ik <= i__2; ik += 2) {
7597 for (ir = 1; ir <= i__3; ++ir) {
7598 bidon += somtab[ir + nd * somtab_dim1] * gsstab[ir + igss *
7602 crvjac[ik + nd * crvjac_dim1] = bidon;
7606 /* --------------- Calculate the coefficients of uneven degree ----------
7610 for (ik = ikdeb + 1; ik <= i__2; ik += 2) {
7614 for (ir = 1; ir <= i__3; ++ir) {
7615 bidon += diftab[ir + nd * diftab_dim1] * gsstab[ir + igss *
7619 crvjac[ik + nd * crvjac_dim1] = bidon;
7626 /* ------- Add terms connected to the supplementary root (0.D0) ------ */
7627 /* ----------- of Legendre polynom of uneven degree NBPNTS -----------
7630 if (*nbpnts % 2 == 0) {
7634 for (nd = 1; nd <= i__1; ++nd) {
7636 for (ik = ikdeb; ik <= i__2; ik += 2) {
7638 crvjac[ik + nd * crvjac_dim1] += somtab[nd * somtab_dim1] *
7639 gsstab[igss * gsstab_dim1];
7645 /* ------------------------------ The end -------------------------------
7650 AdvApp2Var_SysBase::mgsomsg_("MMMAPCO", 7L);
7654 //=======================================================================
7655 //function : mmaperm_
7657 //=======================================================================
7658 int mmaperm_(integer *ncofmx,
7666 /* System generated locals */
7667 integer crvjac_dim1, crvjac_offset, i__1, i__2;
7669 /* Local variables */
7671 integer i__, ia, nd, ncfcut, ibb;
7674 /* **********************************************************************
7679 /* Calculate the square root of the average quadratic error */
7680 /* of approximation done when only the */
7681 /* first NCFNEW coefficients of a curve of degree NCOEFF-1 */
7682 /* written in NORMALIZED Jacobi base of order 2*(IORDRE+1) are preserved. */
7686 /* LEGENDRE,POLYGONE,APPROXIMATION,ERREUR. */
7688 /* INPUT ARGUMENTS : */
7689 /* ------------------ */
7690 /* NCOFMX : Maximum degree of the curve. */
7691 /* NDIM : Dimension of the space. */
7692 /* NCOEFF : Degree +1 of the curve. */
7693 /* IORDRE : Order of constraint of continuity at the extremities. */
7694 /* CRVJAC : The curve the degree which of will be lowered. */
7695 /* NCFNEW : Degree +1 of the resulting polynom. */
7697 /* OUTPUT ARGUMENTS : */
7698 /* ------------------- */
7699 /* ERRMOY : Average precision of approximation. */
7701 /* COMMONS USED : */
7702 /* ---------------- */
7704 /* REFERENCES CALLED : */
7705 /* ----------------------- */
7707 /* DESCRIPTION/NOTES/LIMITATIONS : */
7708 /* ----------------------------------- */
7710 /* ***********************************************************************
7713 /* Name of the routine */
7715 /* Parameter adjustments */
7716 crvjac_dim1 = *ncofmx;
7717 crvjac_offset = crvjac_dim1 + 1;
7718 crvjac -= crvjac_offset;
7721 ibb = AdvApp2Var_SysBase::mnfndeb_();
7723 AdvApp2Var_SysBase::mgenmsg_("MMAPERM", 7L);
7726 /* --------- Minimum degree that can be reached : Stop at 1 or IA -------
7729 ia = (*iordre + 1) << 1;
7731 if (*ncfnew + 1 > ncfcut) {
7732 ncfcut = *ncfnew + 1;
7735 /* -------------- Elimination of coefficients of high degree ------------ */
7736 /* ----------- Loop on the series of Jacobi :NCFCUT --> NCOEFF --------- */
7741 for (nd = 1; nd <= i__1; ++nd) {
7743 for (i__ = ncfcut; i__ <= i__2; ++i__) {
7744 bidj = crvjac[i__ + nd * crvjac_dim1];
7751 /* ----------- Square Root of average quadratic error e -----------
7755 *errmoy = sqrt(bid);
7757 /* ------------------------------- The end ------------------------------
7761 AdvApp2Var_SysBase::mgsomsg_("MMAPERM", 7L);
7765 //=======================================================================
7766 //function : mmapptt_
7768 //=======================================================================
7769 int AdvApp2Var_ApproxF2var::mmapptt_(const integer *ndgjac,
7770 const integer *nbpnts,
7771 const integer *jordre,
7775 /* System generated locals */
7776 integer cgauss_dim1, i__1;
7778 /* Local variables */
7779 integer kjac, iptt, ipdb0, infdg, iptdb, mxjac, ilong, ibb;
7781 /* **********************************************************************
7786 /* Load the elements required for integration by */
7787 /* Gauss method to obtain the coefficients in the base of
7788 /* Legendre of the approximation by the least squares of a */
7789 /* function. The elements are stored in commons MMAPGSS */
7790 /* (case without constraint), MMAPGS0 (constraints C0), MMAPGS1 */
7791 /* (constraints C1) and MMAPGS2 (constraints C2). */
7795 /* INTEGRATION,GAUSS,JACOBI */
7797 /* INPUT ARGUMENTS : */
7798 /* ------------------ */
7799 /* NDGJAC : Max degree of the polynom of approximation. */
7800 /* The representation in orthogonal base goes from degree
7801 /* 0 to degree NDGJAC-2*(JORDRE+1). The polynomial base */
7802 /* is the base of Jacobi of order -1 (Legendre), 0, 1 and 2 */
7803 /* NBPNTS : Degree of the polynom of Legendre on the roots which of */
7804 /* are calculated the coefficients of integration by the */
7805 /* method of Gauss. It is required that NBPNTS=8,10,15,20,25, */
7806 /* 30,40,50 or 61 and NDGJAC < NBPNTS. */
7807 /* JORDRE : Order of the base of Jacobi (-1,0,1 or 2). Corresponds */
7808 /* to step of constraints C0,C1 or C2. */
7810 /* OUTPUT ARGUMENTS : */
7811 /* ------------------- */
7812 /* CGAUSS(i,k) : Table of coefficients of integration by */
7813 /* Gauss method : i varies from 0 to the integer part */
7814 /* of NBPNTS/2 and k varies from 0 to NDGJAC-2*(JORDRE+1). */
7815 /* These are the coeff. of integration associated to */
7816 /* positive roots of the polynom of Legendre of degree */
7817 /* NBPNTS. CGAUSS(0,k) contains coeff. */
7818 /* of integration associated to root t = 0 when */
7819 /* NBPNTS is uneven. */
7820 /* IERCOD : Error code. */
7822 /* = 11 NBPNTS is not 8,10,15,20,25,30,40,50 or 61. */
7823 /* = 21 JORDRE is not -1,0,1 or 2. */
7824 /* = 31 NDGJAC is too great or too small. */
7826 /* COMMONS USED : */
7827 /* ---------------- */
7828 /* MMAPGSS,MMAPGS0,MMAPGS1,MMAPGS2. */
7829 /* ***********************************************************************
7831 /* Parameter adjustments */
7832 cgauss_dim1 = *nbpnts / 2 + 1;
7835 ibb = AdvApp2Var_SysBase::mnfndeb_();
7837 AdvApp2Var_SysBase::mgenmsg_("MMAPPTT", 7L);
7841 /* ------------------- Tests on the validity of inputs ----------------
7844 infdg = (*jordre + 1) << 1;
7845 if (*nbpnts != 8 && *nbpnts != 10 && *nbpnts != 15 && *nbpnts != 20 && *
7846 nbpnts != 25 && *nbpnts != 30 && *nbpnts != 40 && *nbpnts != 50 &&
7851 if (*jordre < -1 || *jordre > 2) {
7855 if (*ndgjac >= *nbpnts || *ndgjac < infdg) {
7859 /* --------------- Calculation of the start pointer following NBPNTS -----------
7864 iptdb += (8 - infdg) << 2;
7867 iptdb += (10 - infdg) * 5;
7870 iptdb += (15 - infdg) * 7;
7873 iptdb += (20 - infdg) * 10;
7876 iptdb += (25 - infdg) * 12;
7879 iptdb += (30 - infdg) * 15;
7882 iptdb += (40 - infdg) * 20;
7885 iptdb += (50 - infdg) * 25;
7890 ipdb0 = ipdb0 + (14 - infdg) / 2 + 1;
7893 ipdb0 = ipdb0 + (24 - infdg) / 2 + 1;
7896 /* ------------------ Choice of the common depending on JORDRE -------------
7899 if (*jordre == -1) {
7912 /* ---------------- Common MMAPGSS (case without constraints) ----------------
7916 ilong = *nbpnts / 2 << 3;
7918 for (kjac = 0; kjac <= i__1; ++kjac) {
7919 iptt = iptdb + kjac * (*nbpnts / 2) + 1;
7920 AdvApp2Var_SysBase::mcrfill_(&ilong,
7921 &mmapgss_.gslxjs[iptt - 1],
7922 &cgauss[kjac * cgauss_dim1 + 1]);
7925 /* --> Case when the number of points is uneven. */
7926 if (*nbpnts % 2 == 1) {
7929 for (kjac = 0; kjac <= i__1; kjac += 2) {
7930 cgauss[kjac * cgauss_dim1] = mmapgss_.gsl0js[iptt - 1];
7935 for (kjac = 1; kjac <= i__1; kjac += 2) {
7936 cgauss[kjac * cgauss_dim1] = 0.;
7942 /* ---------------- Common MMAPGS0 (case with constraints C0) -------------
7946 mxjac = *ndgjac - infdg;
7947 ilong = *nbpnts / 2 << 3;
7949 for (kjac = 0; kjac <= i__1; ++kjac) {
7950 iptt = iptdb + kjac * (*nbpnts / 2) + 1;
7951 AdvApp2Var_SysBase::mcrfill_(&ilong,
7952 &mmapgs0_.gslxj0[iptt - 1],
7953 &cgauss[kjac * cgauss_dim1 + 1]);
7956 /* --> Case when the number of points is uneven. */
7957 if (*nbpnts % 2 == 1) {
7960 for (kjac = 0; kjac <= i__1; kjac += 2) {
7961 cgauss[kjac * cgauss_dim1] = mmapgs0_.gsl0j0[iptt - 1];
7966 for (kjac = 1; kjac <= i__1; kjac += 2) {
7967 cgauss[kjac * cgauss_dim1] = 0.;
7973 /* ---------------- Common MMAPGS1 (case with constraints C1) -------------
7977 mxjac = *ndgjac - infdg;
7978 ilong = *nbpnts / 2 << 3;
7980 for (kjac = 0; kjac <= i__1; ++kjac) {
7981 iptt = iptdb + kjac * (*nbpnts / 2) + 1;
7982 AdvApp2Var_SysBase::mcrfill_(&ilong,
7983 &mmapgs1_.gslxj1[iptt - 1],
7984 &cgauss[kjac * cgauss_dim1 + 1]);
7987 /* --> Case when the number of points is uneven. */
7988 if (*nbpnts % 2 == 1) {
7991 for (kjac = 0; kjac <= i__1; kjac += 2) {
7992 cgauss[kjac * cgauss_dim1] = mmapgs1_.gsl0j1[iptt - 1];
7997 for (kjac = 1; kjac <= i__1; kjac += 2) {
7998 cgauss[kjac * cgauss_dim1] = 0.;
8004 /* ---------------- Common MMAPGS2 (case with constraints C2) -------------
8008 mxjac = *ndgjac - infdg;
8009 ilong = *nbpnts / 2 << 3;
8011 for (kjac = 0; kjac <= i__1; ++kjac) {
8012 iptt = iptdb + kjac * (*nbpnts / 2) + 1;
8013 AdvApp2Var_SysBase::mcrfill_(&ilong,
8014 &mmapgs2_.gslxj2[iptt - 1],
8015 &cgauss[kjac * cgauss_dim1 + 1]);
8018 /* --> Cas of uneven number of points. */
8019 if (*nbpnts % 2 == 1) {
8022 for (kjac = 0; kjac <= i__1; kjac += 2) {
8023 cgauss[kjac * cgauss_dim1] = mmapgs2_.gsl0j2[iptt - 1];
8028 for (kjac = 1; kjac <= i__1; kjac += 2) {
8029 cgauss[kjac * cgauss_dim1] = 0.;
8035 /* ------------------------- Return the error code --------------
8037 /* --> NBPNTS is not OK */
8041 /* --> JORDRE is not OK */
8045 /* --> NDGJAC is not OK */
8050 /* -------------------------------- The end -----------------------------
8055 AdvApp2Var_SysBase::maermsg_("MMAPPTT", iercod, 7L);
8058 AdvApp2Var_SysBase::mgsomsg_("MMAPPTT", 7L);
8064 //=======================================================================
8065 //function : mmjacpt_
8067 //=======================================================================
8068 int mmjacpt_(const integer *ndimen,
8069 const integer *ncoefu,
8070 const integer *ncoefv,
8071 const integer *iordru,
8072 const integer *iordrv,
8073 const doublereal *ptclgd,
8077 /* System generated locals */
8078 integer ptccan_dim1, ptccan_dim2, ptccan_offset, ptclgd_dim1, ptclgd_dim2,
8079 ptclgd_offset, ptcaux_dim1, ptcaux_dim2, ptcaux_dim3,
8080 ptcaux_offset, i__1, i__2, i__3;
8082 /* Local variables */
8083 integer kdim, nd, ii, jj, ibb;
8085 /* ***********************************************************************
8090 /* Passage from canonical to Jacobi base for a */
8091 /* "square" in a space of arbitrary dimension. */
8095 /* SMOOTHING,BASE,LEGENDRE */
8098 /* INPUT ARGUMENTS : */
8099 /* ------------------ */
8100 /* NDIMEN : Dimension of the space. */
8101 /* NCOEFU : Degree+1 by U. */
8102 /* NCOEFV : Degree+1 by V. */
8103 /* IORDRU : Order of Jacobi polynoms by U. */
8104 /* IORDRV : Order of Jacobi polynoms by V. */
8105 /* PTCLGD : The square in the Jacobi base. */
8107 /* OUTPUT ARGUMENTS : */
8108 /* ------------------- */
8109 /* PTCAUX : Auxilliary space. */
8110 /* PTCCAN : The square in the canonic base (-1,1) */
8112 /* COMMONS USED : */
8113 /* ---------------- */
8115 /* APPLIED REFERENCES : */
8116 /* ----------------------- */
8118 /* DESCRIPTION/NOTES/LIMITATIONS : */
8119 /* ----------------------------------- */
8120 /* Cancels and replaces MJACPC */
8122 /* *********************************************************************
8124 /* Name of the routine */
8127 /* Parameter adjustments */
8128 ptccan_dim1 = *ncoefu;
8129 ptccan_dim2 = *ncoefv;
8130 ptccan_offset = ptccan_dim1 * (ptccan_dim2 + 1) + 1;
8131 ptccan -= ptccan_offset;
8132 ptcaux_dim1 = *ncoefv;
8133 ptcaux_dim2 = *ncoefu;
8134 ptcaux_dim3 = *ndimen;
8135 ptcaux_offset = ptcaux_dim1 * (ptcaux_dim2 * (ptcaux_dim3 + 1) + 1) + 1;
8136 ptcaux -= ptcaux_offset;
8137 ptclgd_dim1 = *ncoefu;
8138 ptclgd_dim2 = *ncoefv;
8139 ptclgd_offset = ptclgd_dim1 * (ptclgd_dim2 + 1) + 1;
8140 ptclgd -= ptclgd_offset;
8143 ibb = AdvApp2Var_SysBase::mnfndeb_();
8145 AdvApp2Var_SysBase::mgenmsg_("MMJACPT", 7L);
8148 /* Passage into canonical by u. */
8150 kdim = *ndimen * *ncoefv;
8151 AdvApp2Var_MathBase::mmjaccv_(ncoefu,
8154 &ptclgd[ptclgd_offset],
8155 &ptcaux[ptcaux_offset],
8156 &ptccan[ptccan_offset]);
8158 /* Swapping of u and v. */
8161 for (nd = 1; nd <= i__1; ++nd) {
8163 for (jj = 1; jj <= i__2; ++jj) {
8165 for (ii = 1; ii <= i__3; ++ii) {
8166 ptcaux[jj + (ii + (nd + ptcaux_dim3) * ptcaux_dim2) *
8167 ptcaux_dim1] = ptccan[ii + (jj + nd * ptccan_dim2) *
8176 /* Passage into canonical by v. */
8178 kdim = *ndimen * *ncoefu;
8179 AdvApp2Var_MathBase::mmjaccv_(ncoefv,
8182 &ptcaux[((ptcaux_dim3 + 1) * ptcaux_dim2 + 1) * ptcaux_dim1 + 1],
8183 &ptccan[ptccan_offset],
8184 &ptcaux[(((ptcaux_dim3 << 1) + 1) * ptcaux_dim2 + 1) * ptcaux_dim1 + 1]);
8186 /* Swapping of u and v. */
8189 for (nd = 1; nd <= i__1; ++nd) {
8191 for (jj = 1; jj <= i__2; ++jj) {
8193 for (ii = 1; ii <= i__3; ++ii) {
8194 ptccan[ii + (jj + nd * ptccan_dim2) * ptccan_dim1] = ptcaux[
8195 jj + (ii + (nd + (ptcaux_dim3 << 1)) * ptcaux_dim2) *
8204 /* ---------------------------- THAT'S ALL FOLKS ------------------------
8208 AdvApp2Var_SysBase::mgsomsg_("MMJACPT", 7L);