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,
280 static integer c__1 = 1;
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 static integer nroo2, ncfhe, nd, ii, kk;
290 static integer ibb, kkm, kkp;
291 static doublereal bid1, bid2, bid3;
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 static integer nd, ii, jj, ibb;
509 static doublereal bid;
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 static integer ideb, ifin, nroo2, ideru, iderv;
647 static doublereal renor;
648 static integer ii, nd, ibb, iim, nbp, iip;
649 static doublereal bid1, bid2;
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 static 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 static doublereal rider, riord;
1452 static integer nd, ibb;
1453 static doublereal bid;
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 static doublereal alinu, blinu, alinv, blinv;
1585 static integer ii, ibb;
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 */
1703 static logical ldbg;
1704 static integer ndgu, ndgv;
1705 static doublereal bidu1, bidu2, bidv1, bidv2;
1706 static integer ioru1, iorv1, ii, nd, jj, ku, kv;
1707 static doublereal cnt1, cnt2, cnt3, cnt4;
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 */
1868 static logical ldbg;
1869 static integer ndgv1, ndgv2, ii, jj, nd, kk;
1870 static 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 */
2014 static logical ldbg;
2015 static integer ndgu1, ndgu2, ii, jj, nd, kk;
2016 static 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 */
2160 static logical ldbg;
2161 static 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,
2313 static integer c__1 = 1;
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 static integer ncfhu, ncfhv, nuroo, nvroo, nd, ii, jj, kk, ll, ibb, kkm,
2329 static doublereal bid1, bid2, bid3, bid4;
2330 static 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,
2673 static integer c__1 = 1;
2674 /* System generated locals */
2675 integer sotbv1_dim1, sotbv1_dim2, sotbv1_offset, sotbv2_dim1, sotbv2_dim2,
2676 sotbv2_offset, ditbv1_dim1, ditbv1_dim2, ditbv1_offset,
2677 ditbv2_dim1, ditbv2_dim2, ditbv2_offset, fpntab_dim1,
2678 fpntab_offset, vhermt_dim1, vhermt_offset, sosotb_dim1,
2679 sosotb_dim2, sosotb_offset, diditb_dim1, diditb_dim2,
2680 diditb_offset, soditb_dim1, soditb_dim2, soditb_offset,
2681 disotb_dim1, disotb_dim2, disotb_offset, i__1, i__2, i__3, i__4;
2683 /* Local variables */
2684 static integer ncfhv, nuroo, nvroo, ii, nd, jj, kk, ibb, jjm, jjp;
2685 static 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,
2985 static integer c__1 = 1;
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 static integer ncfhu, nuroo, nvroo, ii, nd, jj, kk, ibb, kkm, kkp;
2998 static 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,
3311 static integer c__8 = 8;
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 */
3322 static integer ilong;
3323 static intptr_t iofwr;
3324 static doublereal wrkar[1];
3325 static integer iszwr;
3326 static integer ibb, ier;
3327 static integer isz1, isz2, isz3, isz4;
3328 static 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;
3476 ibb = AdvApp2Var_SysBase::mnfndeb_();
3478 AdvApp2Var_SysBase::mgenmsg_("MMA2CDI", 7L);
3482 if (*iordru < -1 || *iordru > 2) {
3485 if (*iordrv < -1 || *iordrv > 2) {
3489 /* ------------------------- Set to zero --------------------------------
3492 ilong = (*nbpntu / 2 + 1) * (*nbpntv / 2 + 1) * *ndimen;
3493 AdvApp2Var_SysBase::mvriraz_(&ilong, &sosotb[sosotb_offset]);
3494 AdvApp2Var_SysBase::mvriraz_(&ilong, &diditb[diditb_offset]);
3495 ilong = *nbpntu / 2 * (*nbpntv / 2) * *ndimen;
3496 AdvApp2Var_SysBase::mvriraz_(&ilong, &soditb[soditb_offset]);
3497 AdvApp2Var_SysBase::mvriraz_(&ilong, &disotb[disotb_offset]);
3498 if (*iordru == -1 && *iordrv == -1) {
3504 isz1 = ((*iordru + 1) << 2) * (*iordru + 1);
3505 isz2 = ((*iordrv + 1) << 2) * (*iordrv + 1);
3506 isz3 = ((*iordru + 1) << 1) * *nbpntu;
3507 isz4 = ((*iordrv + 1) << 1) * *nbpntv;
3508 iszwr = isz1 + isz2 + isz3 + isz4;
3509 AdvApp2Var_SysBase::mcrrqst_(&c__8, &iszwr, wrkar, &iofwr, &ier);
3518 if (*iordru >= 0 && *iordru <= 2) {
3520 /* --- Return 2*(IORDRU+1) coeff of 2*(IORDRU+1) polynoms of Hermite
3523 AdvApp2Var_ApproxF2var::mma1her_(iordru, &wrkar[ipt1], iercod);
3528 /* ---- Subract discretizations of polynoms of constraints
3531 mma2cd3_(ndimen, nbpntu, &urootl[1], nbpntv, iordru, &sotbu1[1], &
3532 sotbu2[1], &ditbu1[1], &ditbu2[1], &wrkar[ipt3], &wrkar[ipt1],
3533 &sosotb[sosotb_offset], &soditb[soditb_offset], &disotb[
3534 disotb_offset], &diditb[diditb_offset]);
3537 if (*iordrv >= 0 && *iordrv <= 2) {
3539 /* --- Return 2*(IORDRV+1) coeff of 2*(IORDRV+1) polynoms of Hermite
3542 AdvApp2Var_ApproxF2var::mma1her_(iordrv, &wrkar[ipt2], iercod);
3547 /* ---- Subtract discretisations of polynoms of constraint
3550 mma2cd2_(ndimen, nbpntu, nbpntv, &vrootl[1], iordrv, &sotbv1[1], &
3551 sotbv2[1], &ditbv1[1], &ditbv2[1], &wrkar[ipt4], &wrkar[ipt2],
3552 &sosotb[sosotb_offset], &soditb[soditb_offset], &disotb[
3553 disotb_offset], &diditb[diditb_offset]);
3556 /* --------------- Subtract constraints of corners ----------------
3559 if (*iordru >= 0 && *iordrv >= 0) {
3560 mma2cd1_(ndimen, nbpntu, &urootl[1], nbpntv, &vrootl[1], iordru,
3561 iordrv, &contr1[contr1_offset], &contr2[contr2_offset], &
3562 contr3[contr3_offset], &contr4[contr4_offset], &wrkar[ipt3], &
3563 wrkar[ipt4], &wrkar[ipt1], &wrkar[ipt2], &sosotb[
3564 sosotb_offset], &soditb[soditb_offset], &disotb[disotb_offset]
3565 , &diditb[diditb_offset]);
3569 /* ------------------------------ The End -------------------------------
3571 /* --> IORDRE is not within the autorised diapason. */
3575 /* --> PB of dynamic allocation. */
3582 AdvApp2Var_SysBase::mcrdelt_(&c__8, &iszwr, wrkar, &iofwr, &ier);
3587 AdvApp2Var_SysBase::maermsg_("MMA2CDI", iercod, 7L);
3589 AdvApp2Var_SysBase::mgsomsg_("MMA2CDI", 7L);
3594 //=======================================================================
3595 //function : mma2ce1_
3597 //=======================================================================
3598 int AdvApp2Var_ApproxF2var::mma2ce1_(integer *numdec,
3626 static integer c__8 = 8;
3628 /* System generated locals */
3629 integer sosotb_dim1, sosotb_dim2, sosotb_offset, disotb_dim1, disotb_dim2,
3630 disotb_offset, soditb_dim1, soditb_dim2, soditb_offset,
3631 diditb_dim1, diditb_dim2, diditb_offset, patjac_dim1, patjac_dim2,
3634 /* Local variables */
3635 static logical ldbg;
3636 static intptr_t iofwr;
3637 static doublereal wrkar[1];
3638 static integer iszwr;
3640 static integer isz1, isz2, isz3, isz4, isz5, isz6, isz7;
3641 static intptr_t ipt1, ipt2, ipt3, ipt4, ipt5, ipt6, ipt7;
3645 /* **********************************************************************
3650 /* Calculation of coefficients of polynomial approximation of degree */
3651 /* (NDJACU,NDJACV) of a function F(u,v), starting from its */
3652 /* discretization on roots of Legendre polynom of degree */
3653 /* NBPNTU by U and NBPNTV by V. */
3657 /* TOUS,AB_SPECIFI::FONCTION&,APPROXIMATION,&POLYNOME,&ERREUR */
3659 /* INPUT ARGUMENTS : */
3660 /* ------------------ */
3661 /* NUMDEC: Indicates if it is POSSIBLE to cut function F(u,v). */
3662 /* = 5, It is POSSIBLE to cut by U or by V or in both directions simultaneously. */
3663 /* = 4, It is POSSIBLE to cut by U or by V BUT NOT in both */
3664 /* directions simultaneously (cutting by V is preferable). */
3665 /* = 3, It is POSSIBLE to cut by U or by V BUT NOT in both */
3666 /* directions simultaneously (cutting by U is preferable). */
3667 /* = 2, It is POSSIBLE to cut only by V (i.e. insert parameter */
3668 /* of cutting Vj). */
3669 /* = 1, It is POSSIBLE to cut only by U (i.e. insert parameter */
3670 /* of cutting Ui). */
3671 /* = 0, It is not POSSIBLE to cut anything */
3672 /* NDIMEN: Dimension of the space. */
3673 /* NBSESP: Nb of independent sub-spaces on which the errors are calculated. */
3674 /* NDIMSE: Table of dimensions of each of sub-spaces. */
3675 /* NDMINU: Minimum degree by U to be preserved for the approximation. */
3676 /* NDMINV: Minimum degree by V to be preserved for the approximation. */
3677 /* NDGULI: Limit of nb of coefficients by U of the solution. */
3678 /* NDGVLI: Limit of nb of coefficients by V of the solution. */
3679 /* NDJACU: Max degree of the polynom of approximation by U. */
3680 /* The representation in the orthogonal base starts from degree */
3681 /* 0 to degree NDJACU-2*(IORDRU+1). The polynomial base is the base of */
3682 /* Jacobi of order -1 (Legendre), 0, 1 or 2. */
3683 /* It is required that 2*IORDRU+1 <= NDMINU <= NDGULI < NDJACU */
3684 /* NDJACV: Max degree of the polynom of approximation by V. */
3685 /* The representation in the orthogonal base starts from degree */
3686 /* 0 to degree NDJACV-2*(IORDRV+1). The polynomial base is */
3687 /* the base of Jacobi of order -1 (Legendre), 0, 1 or 2 */
3688 /* It is required that 2*IORDRV+1 <= NDMINV <= NDGVLI < NDJACV */
3689 /* IORDRU: Order of the Jacobi base (-1,0,1 or 2) by U. Corresponds */
3690 /* to the step of constraints C0, C1 or C2. */
3691 /* IORDRV: Order of the Jacobi base (-1,0,1 or 2) by U. Corresponds */
3692 /* to the step of constraints C0, C1 or C2. */
3693 /* NBPNTU: Degree of Legendre polynom on the roots which of are */
3694 /* calculated the coefficients of integration by u */
3695 /* by Gauss method. It is required that NBPNTU = 30, 40, */
3696 /* 50 or 61 and NDJACU-2*(IORDRU+1) < NBPNTU. */
3697 /* NBPNTV: Degree of Legendre polynom on the roots which of are */
3698 /* calculated the coefficients of integration by u */
3699 /* by Gauss method. It is required that NBPNTV = 30, 40, */
3700 /* 50 or 61 and NDJACV-2*(IORDRV+1) < NBPNTV. */
3701 /* EPSAPR: Table of NBSESP tolerances imposed on each sub-spaces. */
3702 /* SOSOTB: Table of F(ui,vj) + F(ui,-vj) + F(-ui,vj) + F(-ui,-vj) */
3703 /* with ui and vj - positive roots of the Legendre polynom */
3704 /* of degree NBPNTU and NBPNTV respectively. Additionally, */
3705 /* table SOSOTB(0,j) contains F(0,vj) + F(0,-vj), */
3706 /* table SOSOTB(i,0) contains F(ui,0) + F(-ui,0) and */
3707 /* SOSOTB(0,0) contains F(0,0). */
3708 /* DISOTB: Table of F(ui,vj) + F(ui,-vj) - F(-ui,vj) - F(-ui,-vj) */
3709 /* with ui and vj positive roots of Legendre polynom */
3710 /* of degree NBPNTU and NBPNTV respectively. */
3711 /* SODITB: Table of F(ui,vj) - F(ui,-vj) + F(-ui,vj) - F(-ui,-vj) */
3712 /* with ui and vj positive roots of Legendre polynom */
3713 /* of degree NBPNTU and NBPNTV respectively. */
3714 /* DIDITB: Table of F(ui,vj) - F(ui,-vj) - F(-ui,vj) + F(-ui,-vj) */
3715 /* with ui and vj positive roots of Legendre polynom */
3716 /* of degree NBPNTU and NBPNTV respectively. Additionally, */
3717 /* table DIDITB(0,j) contains F(0,vj) - F(0,-vj), */
3718 /* and table DIDITB(i,0) contains F(ui,0) - F(-ui,0). */
3720 /* OUTPUT ARGUMENTS */
3721 /* --------------- */
3722 /* PATJAC: Table of coefficients of polynom P(u,v) of approximation */
3723 /* of F(u,v) with eventually taking into account of */
3724 /* constraints. P(u,v) is of degree (NDJACU,NDJACV). */
3725 /* This table contains other coeff if ITYDEC = 0. */
3726 /* ERRMAX: For 1<=i<=NBSESP, ERRMAX(i) contains max errors */
3727 /* on each of sub-spaces SI ITYDEC = 0. */
3728 /* ERRMOY: Contains average errors for each of NBSESP sub-spaces SI ITYDEC = 0. */
3729 /* NDEGPU: Degree by U for square PATJAC. Valable if ITYDEC=0. */
3730 /* NDEGPV: Degree by V for square PATJAC. Valable if ITYDEC=0. */
3731 /* ITYDEC: Shows if it is NECESSARY to cut again function F(u,v). */
3732 /* = 0, it is not NECESSARY to cut anything, PATJAC is OK. */
3733 /* = 1, it is NECESSARY to cut only by U (i.e. insert parameter of cutting Ui). */
3734 /* = 2, it is NECESSARY to cut only by V (i.e. insert parameter of cutting Vj). */
3735 /* = 3, it is NECESSARY to cut both by U AND by V. */
3736 /* IERCOD: Error code. */
3737 /* = 0, Everything is OK. */
3738 /* = -1, There is the best possible solution, but the */
3739 /* user tolerance is not satisfactory (3*only) */
3740 /* = 1, Incoherent entries. */
3742 /* COMMONS USED : */
3743 /* ---------------- */
3745 /* REFERENCES CALLED : */
3746 /* --------------------- */
3748 /* DESCRIPTION/NOTES/LIMITATIONS : */
3749 /* ------------------------------- */
3752 /* **********************************************************************
3754 /* Name of the routine */
3757 /* --------------------------- Initialisations --------------------------
3760 /* Parameter adjustments */
3765 patjac_dim1 = *ndjacu + 1;
3766 patjac_dim2 = *ndjacv + 1;
3767 patjac_offset = patjac_dim1 * patjac_dim2;
3768 patjac -= patjac_offset;
3769 diditb_dim1 = *nbpntu / 2 + 1;
3770 diditb_dim2 = *nbpntv / 2 + 1;
3771 diditb_offset = diditb_dim1 * diditb_dim2;
3772 diditb -= diditb_offset;
3773 soditb_dim1 = *nbpntu / 2;
3774 soditb_dim2 = *nbpntv / 2;
3775 soditb_offset = soditb_dim1 * (soditb_dim2 + 1) + 1;
3776 soditb -= soditb_offset;
3777 disotb_dim1 = *nbpntu / 2;
3778 disotb_dim2 = *nbpntv / 2;
3779 disotb_offset = disotb_dim1 * (disotb_dim2 + 1) + 1;
3780 disotb -= disotb_offset;
3781 sosotb_dim1 = *nbpntu / 2 + 1;
3782 sosotb_dim2 = *nbpntv / 2 + 1;
3783 sosotb_offset = sosotb_dim1 * sosotb_dim2;
3784 sosotb -= sosotb_offset;
3787 ldbg = AdvApp2Var_SysBase::mnfndeb_() >= 3;
3789 AdvApp2Var_SysBase::mgenmsg_("MMA2CE1", 7L);
3794 isz1 = (*nbpntu / 2 + 1) * (*ndjacu - ((*iordru + 1) << 1) + 1);
3795 isz2 = (*nbpntv / 2 + 1) * (*ndjacv - ((*iordrv + 1) << 1) + 1);
3796 isz3 = (*nbpntv / 2 + 1) * (*ndjacu - ((*iordru + 1) << 1) + 1) * *ndimen;
3797 isz4 = *nbpntv / 2 * (*ndjacu - ((*iordru + 1) << 1) + 1) * *ndimen;
3798 isz5 = *ndjacu + 1 - ((*iordru + 1) << 1);
3799 isz6 = *ndjacv + 1 - ((*iordrv + 1) << 1);
3800 isz7 = *ndimen << 2;
3801 iszwr = isz1 + isz2 + isz3 + isz4 + isz5 + isz6 + isz7;
3802 AdvApp2Var_SysBase::mcrrqst_(&c__8, &iszwr, wrkar, &iofwr, &ier);
3814 /* ----------------- Return Gauss coefficients of integration ----------------
3817 AdvApp2Var_ApproxF2var::mmapptt_(ndjacu, nbpntu, iordru, &wrkar[ipt1], iercod);
3821 AdvApp2Var_ApproxF2var::mmapptt_(ndjacv, nbpntv, iordrv, &wrkar[ipt2], iercod);
3826 /* ------------------- Return max polynoms of Jacobi ------------
3829 AdvApp2Var_ApproxF2var::mma2jmx_(ndjacu, iordru, &wrkar[ipt5]);
3830 AdvApp2Var_ApproxF2var::mma2jmx_(ndjacv, iordrv, &wrkar[ipt6]);
3832 /* ------ Calculate the coefficients and their contribution to the error ----
3835 mma2ce2_(numdec, ndimen, nbsesp, &ndimse[1], ndminu, ndminv, ndguli,
3836 ndgvli, ndjacu, ndjacv, iordru, iordrv, nbpntu, nbpntv, &epsapr[1]
3837 , &sosotb[sosotb_offset], &disotb[disotb_offset], &soditb[
3838 soditb_offset], &diditb[diditb_offset], &wrkar[ipt1], &wrkar[ipt2]
3839 , &wrkar[ipt5], &wrkar[ipt6], &wrkar[ipt7], &wrkar[ipt3], &wrkar[
3840 ipt4], &patjac[patjac_offset], &errmax[1], &errmoy[1], ndegpu,
3841 ndegpv, itydec, iercod);
3847 /* ------------------------------ The end -------------------------------
3856 AdvApp2Var_SysBase::mcrdelt_(&c__8, &iszwr, wrkar, &iofwr, &ier);
3861 AdvApp2Var_SysBase::maermsg_("MMA2CE1", iercod, 7L);
3863 AdvApp2Var_SysBase::mgsomsg_("MMA2CE1", 7L);
3868 //=======================================================================
3869 //function : mma2ce2_
3871 //=======================================================================
3872 int mma2ce2_(integer *numdec,
3907 /* System generated locals */
3908 integer sosotb_dim1, sosotb_dim2, sosotb_offset, disotb_dim1, disotb_dim2,
3909 disotb_offset, soditb_dim1, soditb_dim2, soditb_offset,
3910 diditb_dim1, diditb_dim2, diditb_offset, gssutb_dim1, gssvtb_dim1,
3911 chpair_dim1, chpair_dim2, chpair_offset, chimpr_dim1,
3912 chimpr_dim2, chimpr_offset, patjac_dim1, patjac_dim2,
3913 patjac_offset, vecerr_dim1, vecerr_offset, i__1, i__2, i__3, i__4;
3915 /* Local variables */
3916 static logical ldbg;
3917 static integer idim, igsu, minu, minv, maxu, maxv, igsv;
3918 static doublereal vaux[3];
3919 static integer i2rdu, i2rdv, ndses, nd, ii, jj, kk, nu, nv;
3920 static doublereal zu, zv;
3921 static integer nu1, nv1;
3923 /* **********************************************************************
3927 /* Calculation of coefficients of polynomial approximation of degree */
3928 /* (NDJACU,NDJACV) of a function F(u,v), starting from its */
3929 /* discretization on roots of Legendre polynom of degree */
3930 /* NBPNTU by U and NBPNTV by V. */
3934 /* TOUS,AB_SPECIFI::FONCTION&,APPROXIMATION,&COEFFICIENT,&POLYNOME */
3936 /* INPUT ARGUMENTS : */
3937 /* ------------------ */
3938 /* NUMDEC: Indicates if it is POSSIBLE to cut function F(u,v). */
3939 /* = 5, It is POSSIBLE to cut by U or by V or in both directions simultaneously. */
3940 /* = 4, It is POSSIBLE to cut by U or by V BUT NOT in both */
3941 /* directions simultaneously (cutting by V is preferable). */
3942 /* = 3, It is POSSIBLE to cut by U or by V BUT NOT in both */
3943 /* directions simultaneously (cutting by U is preferable). */
3944 /* = 2, It is POSSIBLE to cut only by V (i.e. insert parameter */
3945 /* of cutting Vj). */
3946 /* = 1, It is POSSIBLE to cut only by U (i.e. insert parameter */
3947 /* of cutting Ui). */
3948 /* = 0, It is not POSSIBLE to cut anything */
3949 /* NDIMEN: Total dimension of the space. */
3950 /* NBSESP: Nb of independent sub-spaces on which the errors are calculated. */
3951 /* NDIMSE: Table of dimensions of each of sub-spaces. */
3952 /* NDMINU: Minimum degree by U to be preserved for the approximation. */
3953 /* NDMINV: Minimum degree by V to be preserved for the approximation. */
3954 /* NDGULI: Limit of nb of coefficients by U of the solution. */
3955 /* NDGVLI: Limit of nb of coefficients by V of the solution. */
3956 /* NDJACU: Max degree of the polynom of approximation by U. */
3957 /* The representation in the orthogonal base starts from degree */
3958 /* 0 to degree NDJACU-2*(IORDRU+1). The polynomial base is the base of */
3959 /* Jacobi of order -1 (Legendre), 0, 1 or 2. */
3960 /* It is required that 2*IORDRU+1 <= NDMINU <= NDGULI < NDJACU */
3961 /* NDJACV: Max degree of the polynom of approximation by V. */
3962 /* The representation in the orthogonal base starts from degree */
3963 /* 0 to degree NDJACV-2*(IORDRV+1). The polynomial base is */
3964 /* the base of Jacobi of order -1 (Legendre), 0, 1 or 2 */
3965 /* It is required that 2*IORDRV+1 <= NDMINV <= NDGVLI < NDJACV */
3966 /* IORDRU: Order of the Jacobi base (-1,0,1 or 2) by U. Corresponds */
3967 /* to the step of constraints C0, C1 or C2. */
3968 /* IORDRV: Order of the Jacobi base (-1,0,1 or 2) by U. Corresponds */
3969 /* to the step of constraints C0, C1 or C2. */
3970 /* NBPNTU: Degree of Legendre polynom on the roots which of are */
3971 /* calculated the coefficients of integration by u */
3972 /* by Gauss method. It is required that NBPNTU = 30, 40, */
3973 /* 50 or 61 and NDJACU-2*(IORDRU+1) < NBPNTU. */
3974 /* NBPNTV: Degree of Legendre polynom on the roots which of are */
3975 /* calculated the coefficients of integration by u */
3976 /* by Gauss method. It is required that NBPNTV = 30, 40, */
3977 /* 50 or 61 and NDJACV-2*(IORDRV+1) < NBPNTV. */
3978 /* EPSAPR: Table of NBSESP tolerances imposed on each sub-spaces. */
3979 /* SOSOTB: Table of F(ui,vj) + F(ui,-vj) + F(-ui,vj) + F(-ui,-vj) */
3980 /* with ui and vj - positive roots of the Legendre polynom */
3981 /* of degree NBPNTU and NBPNTV respectively. Additionally, */
3982 /* table SOSOTB(0,j) contains F(0,vj) + F(0,-vj), */
3983 /* table SOSOTB(i,0) contains F(ui,0) + F(-ui,0) and */
3984 /* SOSOTB(0,0) contains F(0,0). */
3985 /* DISOTB: Table of F(ui,vj) + F(ui,-vj) - F(-ui,vj) - F(-ui,-vj) */
3986 /* with ui and vj positive roots of Legendre polynom */
3987 /* of degree NBPNTU and NBPNTV respectively. */
3988 /* SODITB: Table of F(ui,vj) - F(ui,-vj) + F(-ui,vj) - F(-ui,-vj) */
3989 /* with ui and vj positive roots of Legendre polynom */
3990 /* of degree NBPNTU and NBPNTV respectively. */
3991 /* DIDITB: Table of F(ui,vj) - F(ui,-vj) - F(-ui,vj) + F(-ui,-vj) */
3992 /* with ui and vj positive roots of Legendre polynom */
3993 /* of degree NBPNTU and NBPNTV respectively. Additionally, */
3994 /* table DIDITB(0,j) contains F(0,vj) - F(0,-vj), */
3995 /* and table DIDITB(i,0) contains F(ui,0) - F(-ui,0). */
3996 /* GSSUTB: Table of coefficients of integration by Gauss method */
3997 /* by U: i varies from 0 to NBPNTU/2 and k varies from 0 to */
3998 /* NDJACU-2*(IORDRU+1). */
3999 /* GSSVTB: Table of coefficients of integration by Gauss method */
4000 /* by V: i varies from 0 to NBPNTV/2 and k varies from 0 to */
4001 /* NDJACV-2*(IORDRV+1). */
4002 /* XMAXJU: Maximum value of Jacobi polynoms of order IORDRU, */
4003 /* from degree 0 to degree NDJACU - 2*(IORDRU+1) */
4004 /* XMAXJV: Maximum value of Jacobi polynoms of order IORDRV, */
4005 /* from degree 0 to degree NDJACV - 2*(IORDRV+1) */
4007 /* OUTPUT ARGUMENTS : */
4008 /* ------------------- */
4009 /* VECERR: Auxiliary table. */
4010 /* CHPAIR: Auxiliary table of terms connected to degree NDJACU by U */
4011 /* to calculate the coeff. of approximation of EVEN degree by V. */
4012 /* CHIMPR: Auxiliary table of terms connected to degree NDJACU by U */
4013 /* to calculate the coeff. of approximation of UNEVEN degree by V. */
4014 /* PATJAC: Table of coefficients of polynom P(u,v) of approximation */
4015 /* of F(u,v) with eventually taking into account of */
4016 /* constraints. P(u,v) is of degree (NDJACU,NDJACV). */
4017 /* This table contains other coeff if ITYDEC = 0. */
4018 /* ERRMAX: For 1<=i<=NBSESP, ERRMAX(i) contains max errors */
4019 /* on each of sub-spaces SI ITYDEC = 0. */
4020 /* ERRMOY: Contains average errors for each of NBSESP sub-spaces SI ITYDEC = 0. */
4021 /* NDEGPU: Degree by U for square PATJAC. Valable if ITYDEC=0. */
4022 /* NDEGPV: Degree by V for square PATJAC. Valable if ITYDEC=0. */
4023 /* ITYDEC: Shows if it is NECESSARY to cut again function F(u,v). */
4024 /* = 0, it is not NECESSARY to cut anything, PATJAC is OK. */
4025 /* = 1, it is NECESSARY to cut only by U (i.e. insert parameter of cutting Ui). */
4026 /* = 2, it is NECESSARY to cut only by V (i.e. insert parameter of cutting Vj). */
4027 /* = 3, it is NECESSARY to cut both by U AND by V. */
4028 /* IERCOD: Error code. */
4029 /* = 0, Everything is OK. */
4030 /* = -1, There is the best possible solution, but the */
4031 /* user tolerance is not satisfactory (3*only) */
4032 /* = 1, Incoherent entries. */
4034 /* COMMONS USED : */
4035 /* ---------------- */
4037 /* REFERENCES CALLED : */
4038 /* --------------------- */
4040 /* DESCRIPTION/NOTES/LIMITATIONS : */
4042 /* **********************************************************************
4044 /* Name of the routine */
4047 /* --------------------------- Initialisations --------------------------
4050 /* Parameter adjustments */
4051 vecerr_dim1 = *ndimen;
4052 vecerr_offset = vecerr_dim1 + 1;
4053 vecerr -= vecerr_offset;
4058 patjac_dim1 = *ndjacu + 1;
4059 patjac_dim2 = *ndjacv + 1;
4060 patjac_offset = patjac_dim1 * patjac_dim2;
4061 patjac -= patjac_offset;
4062 gssutb_dim1 = *nbpntu / 2 + 1;
4063 chimpr_dim1 = *nbpntv / 2;
4064 chimpr_dim2 = *ndjacu - ((*iordru + 1) << 1) + 1;
4065 chimpr_offset = chimpr_dim1 * chimpr_dim2 + 1;
4066 chimpr -= chimpr_offset;
4067 chpair_dim1 = *nbpntv / 2 + 1;
4068 chpair_dim2 = *ndjacu - ((*iordru + 1) << 1) + 1;
4069 chpair_offset = chpair_dim1 * chpair_dim2;
4070 chpair -= chpair_offset;
4071 gssvtb_dim1 = *nbpntv / 2 + 1;
4072 diditb_dim1 = *nbpntu / 2 + 1;
4073 diditb_dim2 = *nbpntv / 2 + 1;
4074 diditb_offset = diditb_dim1 * diditb_dim2;
4075 diditb -= diditb_offset;
4076 soditb_dim1 = *nbpntu / 2;
4077 soditb_dim2 = *nbpntv / 2;
4078 soditb_offset = soditb_dim1 * (soditb_dim2 + 1) + 1;
4079 soditb -= soditb_offset;
4080 disotb_dim1 = *nbpntu / 2;
4081 disotb_dim2 = *nbpntv / 2;
4082 disotb_offset = disotb_dim1 * (disotb_dim2 + 1) + 1;
4083 disotb -= disotb_offset;
4084 sosotb_dim1 = *nbpntu / 2 + 1;
4085 sosotb_dim2 = *nbpntv / 2 + 1;
4086 sosotb_offset = sosotb_dim1 * sosotb_dim2;
4087 sosotb -= sosotb_offset;
4090 ldbg = AdvApp2Var_SysBase::mnfndeb_() >= 3;
4092 AdvApp2Var_SysBase::mgenmsg_("MMA2CE2", 7L);
4094 /* --> A priori everything is OK */
4096 /* --> test of inputs */
4097 if (*numdec < 0 || *numdec > 5) {
4100 if ((*iordru << 1) + 1 > *ndminu) {
4103 if (*ndminu > *ndguli) {
4106 if (*ndguli >= *ndjacu) {
4109 if ((*iordrv << 1) + 1 > *ndminv) {
4112 if (*ndminv > *ndgvli) {
4115 if (*ndgvli >= *ndjacv) {
4118 /* --> A priori, no cuts to be done */
4120 /* --> Min. degrees to return: NDMINU,NDMINV */
4123 /* --> For the moment, max errors are null */
4124 AdvApp2Var_SysBase::mvriraz_(nbsesp, &errmax[1]);
4126 AdvApp2Var_SysBase::mvriraz_(&nd, &vecerr[vecerr_offset]);
4127 /* --> and the square, too. */
4128 nd = (*ndjacu + 1) * (*ndjacv + 1) * *ndimen;
4129 AdvApp2Var_SysBase::mvriraz_(&nd, &patjac[patjac_offset]);
4131 i2rdu = (*iordru + 1) << 1;
4132 i2rdv = (*iordrv + 1) << 1;
4134 /* **********************************************************************
4136 /* -------------------- HERE IT IS POSSIBLE TO CUT ----------------------
4138 /* **********************************************************************
4141 if (*numdec > 0 && *numdec <= 5) {
4143 /* ******************************************************************
4145 /* ---------------------- Calculate coeff of zone 4 -------------
4159 /* ---------------- Calculate the terms connected to degree by U ---------
4163 for (nd = 1; nd <= i__1; ++nd) {
4165 for (kk = minu; kk <= i__2; ++kk) {
4167 mma2cfu_(&kk, nbpntu, nbpntv, &sosotb[nd * sosotb_dim2 *
4168 sosotb_dim1], &disotb[(nd * disotb_dim2 + 1) *
4169 disotb_dim1 + 1], &soditb[(nd * soditb_dim2 + 1) *
4170 soditb_dim1 + 1], &diditb[nd * diditb_dim2 *
4171 diditb_dim1], &gssutb[igsu * gssutb_dim1], &chpair[(
4172 igsu + nd * chpair_dim2) * chpair_dim1], &chimpr[(
4173 igsu + nd * chimpr_dim2) * chimpr_dim1 + 1]);
4179 /* ------------------- Calculate the coefficients of PATJAC ------------
4182 igsu = minu - i2rdu;
4184 for (jj = minv; jj <= i__1; ++jj) {
4187 for (nd = 1; nd <= i__2; ++nd) {
4188 mma2cfv_(&jj, &minu, &maxu, nbpntv, &gssvtb[igsv *
4189 gssvtb_dim1], &chpair[(igsu + nd * chpair_dim2) *
4190 chpair_dim1], &chimpr[(igsu + nd * chimpr_dim2) *
4191 chimpr_dim1 + 1], &patjac[minu + (jj + nd *
4192 patjac_dim2) * patjac_dim1]);
4196 /* ----- Contribution of calculated terms to the approximation error */
4197 /* for terms (I,J) with MINU <= I <= MAXU, J fixe. */
4201 for (nd = 1; nd <= i__2; ++nd) {
4203 mma2er1_(ndjacu, ndjacv, &ndses, &minu, &maxu, &jj, &jj,
4204 iordru, iordrv, xmaxju, xmaxjv, &patjac[idim *
4205 patjac_dim2 * patjac_dim1], &vecerr[vecerr_dim1 + 1],
4206 &vecerr[nd + (vecerr_dim1 << 2)]);
4207 if (vecerr[nd + (vecerr_dim1 << 2)] > epsapr[nd]) {
4216 /* ******************************************************************
4218 /* ---------------------- Calculate the coeff of zone 2 -------------
4221 minu = (*iordru + 1) << 1;
4226 /* --> If zone 2 is empty, pass to zone 3. */
4227 /* VECERR(ND,2) was already set to zero. */
4232 /* ---------------- Calculate the terms connected to degree by U ------------
4236 for (nd = 1; nd <= i__1; ++nd) {
4238 for (kk = minu; kk <= i__2; ++kk) {
4240 mma2cfu_(&kk, nbpntu, nbpntv, &sosotb[nd * sosotb_dim2 *
4241 sosotb_dim1], &disotb[(nd * disotb_dim2 + 1) *
4242 disotb_dim1 + 1], &soditb[(nd * soditb_dim2 + 1) *
4243 soditb_dim1 + 1], &diditb[nd * diditb_dim2 *
4244 diditb_dim1], &gssutb[igsu * gssutb_dim1], &chpair[(
4245 igsu + nd * chpair_dim2) * chpair_dim1], &chimpr[(
4246 igsu + nd * chimpr_dim2) * chimpr_dim1 + 1]);
4252 /* ------------------- Calculate the coefficients of PATJAC ------------
4255 igsu = minu - i2rdu;
4257 for (jj = minv; jj <= i__1; ++jj) {
4260 for (nd = 1; nd <= i__2; ++nd) {
4261 mma2cfv_(&jj, &minu, &maxu, nbpntv, &gssvtb[igsv *
4262 gssvtb_dim1], &chpair[(igsu + nd * chpair_dim2) *
4263 chpair_dim1], &chimpr[(igsu + nd * chimpr_dim2) *
4264 chimpr_dim1 + 1], &patjac[minu + (jj + nd *
4265 patjac_dim2) * patjac_dim1]);
4271 /* -----Contribution of calculated terms to the approximation error */
4272 /* for terms (I,J) with MINU <= I <= MAXU, MINV <= J <= MAXV */
4276 for (nd = 1; nd <= i__1; ++nd) {
4278 mma2er1_(ndjacu, ndjacv, &ndses, &minu, &maxu, &minv, &maxv,
4279 iordru, iordrv, xmaxju, xmaxjv, &patjac[idim *
4280 patjac_dim2 * patjac_dim1], &vecerr[vecerr_dim1 + 1], &
4281 vecerr[nd + (vecerr_dim1 << 1)]);
4286 /* ******************************************************************
4288 /* ---------------------- Calculation of coeff of zone 3 -------------
4294 minv = (*iordrv + 1) << 1;
4297 /* -> If zone 3 is empty, pass to the test of cutting. */
4298 /* VECERR(ND,3) was already set to zero */
4303 /* ----------- The terms connected to the degree by U are already calculated -----
4305 /* ------------------- Calculation of coefficients of PATJAC ------------
4308 igsu = minu - i2rdu;
4310 for (jj = minv; jj <= i__1; ++jj) {
4313 for (nd = 1; nd <= i__2; ++nd) {
4314 mma2cfv_(&jj, &minu, &maxu, nbpntv, &gssvtb[igsv *
4315 gssvtb_dim1], &chpair[(igsu + nd * chpair_dim2) *
4316 chpair_dim1], &chimpr[(igsu + nd * chimpr_dim2) *
4317 chimpr_dim1 + 1], &patjac[minu + (jj + nd *
4318 patjac_dim2) * patjac_dim1]);
4324 /* ----- Contribution of calculated terms to the approximation error
4325 /* for terms (I,J) with MINU <= I <= MAXU, MINV <= J <= MAXV. */
4329 for (nd = 1; nd <= i__1; ++nd) {
4331 mma2er1_(ndjacu, ndjacv, &ndses, &minu, &maxu, &minv, &maxv,
4332 iordru, iordrv, xmaxju, xmaxjv, &patjac[idim *
4333 patjac_dim2 * patjac_dim1], &vecerr[vecerr_dim1 + 1], &
4334 vecerr[nd + vecerr_dim1 * 3]);
4339 /* ******************************************************************
4341 /* --------------------------- Tests of cutting ---------------------
4346 for (nd = 1; nd <= i__1; ++nd) {
4347 vaux[0] = vecerr[nd + (vecerr_dim1 << 1)];
4348 vaux[1] = vecerr[nd + (vecerr_dim1 << 2)];
4349 vaux[2] = vecerr[nd + vecerr_dim1 * 3];
4351 errmax[nd] = AdvApp2Var_MathBase::mzsnorm_(&ii, vaux);
4352 if (errmax[nd] > epsapr[nd]) {
4354 zv = AdvApp2Var_MathBase::mzsnorm_(&ii, vaux);
4355 zu = AdvApp2Var_MathBase::mzsnorm_(&ii, &vaux[1]);
4356 if (zu > epsapr[nd] && zv > epsapr[nd]) {
4368 /* ******************************************************************
4370 /* --- OK, the square is valid, the coeff of zone 1 are calculated
4373 minu = (*iordru + 1) << 1;
4375 minv = (*iordrv + 1) << 1;
4378 /* --> If zone 1 is empty, pass to the calculation of Max and Average error. */
4379 if (minu > maxu || minv > maxv) {
4383 /* ----------- The terms connected to degree by U are already calculated -----
4385 /* ------------------- Calculate the coefficients of PATJAC ------------
4388 igsu = minu - i2rdu;
4390 for (jj = minv; jj <= i__1; ++jj) {
4393 for (nd = 1; nd <= i__2; ++nd) {
4394 mma2cfv_(&jj, &minu, &maxu, nbpntv, &gssvtb[igsv *
4395 gssvtb_dim1], &chpair[(igsu + nd * chpair_dim2) *
4396 chpair_dim1], &chimpr[(igsu + nd * chimpr_dim2) *
4397 chimpr_dim1 + 1], &patjac[minu + (jj + nd *
4398 patjac_dim2) * patjac_dim1]);
4404 /* --------------- Now the degree is maximally lowered --------
4409 i__1 = 1, i__2 = (*iordru << 1) + 1, i__1 = advapp_max(i__1,i__2);
4410 minu = advapp_max(i__1,*ndminu);
4413 i__1 = 1, i__2 = (*iordrv << 1) + 1, i__1 = advapp_max(i__1,i__2);
4414 minv = advapp_max(i__1,*ndminv);
4418 for (nd = 1; nd <= i__1; ++nd) {
4420 if (maxu >= (*iordru + 1) << 1 && maxv >= (*iordrv + 1) << 1) {
4421 mma2er2_(ndjacu, ndjacv, &ndses, &minu, &maxu, &minv, &maxv,
4422 iordru, iordrv, xmaxju, xmaxjv, &patjac[idim *
4423 patjac_dim2 * patjac_dim1], &epsapr[nd], &vecerr[
4424 vecerr_dim1 + 1], &errmax[nd], &nu, &nv);
4432 /* --> Calculate the average error. */
4433 mma2moy_(ndjacu, ndjacv, &ndses, &nu1, ndjacu, &nv1, ndjacv,
4434 iordru, iordrv, &patjac[idim * patjac_dim2 * patjac_dim1],
4437 /* --> Set to 0.D0 the rejected coeffs. */
4438 i__2 = idim + ndses - 1;
4439 for (ii = idim; ii <= i__2; ++ii) {
4441 for (jj = nv1; jj <= i__3; ++jj) {
4443 for (kk = nu1; kk <= i__4; ++kk) {
4444 patjac[kk + (jj + ii * patjac_dim2) * patjac_dim1] =
4453 /* --> Return the nb of coeffs of approximation. */
4454 *ndegpu = advapp_max(*ndegpu,nu);
4455 *ndegpv = advapp_max(*ndegpv,nv);
4460 /* ******************************************************************
4462 /* -------------------- IT IS NOT POSSIBLE TO CUT -------------------
4464 /* ******************************************************************
4468 minu = (*iordru + 1) << 1;
4470 minv = (*iordrv + 1) << 1;
4473 /* ---------------- Calculate the terms connected to the degree by U ------------
4477 for (nd = 1; nd <= i__1; ++nd) {
4479 for (kk = minu; kk <= i__2; ++kk) {
4481 mma2cfu_(&kk, nbpntu, nbpntv, &sosotb[nd * sosotb_dim2 *
4482 sosotb_dim1], &disotb[(nd * disotb_dim2 + 1) *
4483 disotb_dim1 + 1], &soditb[(nd * soditb_dim2 + 1) *
4484 soditb_dim1 + 1], &diditb[nd * diditb_dim2 *
4485 diditb_dim1], &gssutb[igsu * gssutb_dim1], &chpair[(
4486 igsu + nd * chpair_dim2) * chpair_dim1], &chimpr[(
4487 igsu + nd * chimpr_dim2) * chimpr_dim1 + 1]);
4491 /* ---------------------- Calculate all coefficients -------
4494 igsu = minu - i2rdu;
4496 for (jj = minv; jj <= i__2; ++jj) {
4498 mma2cfv_(&jj, &minu, &maxu, nbpntv, &gssvtb[igsv *
4499 gssvtb_dim1], &chpair[(igsu + nd * chpair_dim2) *
4500 chpair_dim1], &chimpr[(igsu + nd * chimpr_dim2) *
4501 chimpr_dim1 + 1], &patjac[minu + (jj + nd *
4502 patjac_dim2) * patjac_dim1]);
4508 /* ----- Contribution of calculated terms to the approximation error
4509 /* for terms (I,J) with MINU <= I <= MAXU, MINV <= J <= MAXV */
4513 for (nd = 1; nd <= i__1; ++nd) {
4515 minu = (*iordru + 1) << 1;
4519 mma2er1_(ndjacu, ndjacv, &ndses, &minu, &maxu, &minv, &maxv,
4520 iordru, iordrv, xmaxju, xmaxjv, &patjac[idim *
4521 patjac_dim2 * patjac_dim1], &vecerr[vecerr_dim1 + 1], &
4525 minv = (*iordrv + 1) << 1;
4528 mma2er1_(ndjacu, ndjacv, &ndses, &minu, &maxu, &minv, &maxv,
4529 iordru, iordrv, xmaxju, xmaxjv, &patjac[idim *
4530 patjac_dim2 * patjac_dim1], &vecerr[vecerr_dim1 + 1],
4534 /* ---------------------------- IF ERRMAX > EPSAPR, stop --------
4537 if (errmax[nd] > epsapr[nd]) {
4542 /* ------------- Otherwise, try to remove again the coeff
4547 i__2 = 1, i__3 = (*iordru << 1) + 1, i__2 = advapp_max(i__2,i__3);
4548 minu = advapp_max(i__2,*ndminu);
4551 i__2 = 1, i__3 = (*iordrv << 1) + 1, i__2 = advapp_max(i__2,i__3);
4552 minv = advapp_max(i__2,*ndminv);
4554 if (maxu >= (*iordru + 1) << 1 && maxv >= (*iordrv + 1) << 1) {
4555 mma2er2_(ndjacu, ndjacv, &ndses, &minu, &maxu, &minv, &
4556 maxv, iordru, iordrv, xmaxju, xmaxjv, &patjac[
4557 idim * patjac_dim2 * patjac_dim1], &epsapr[nd], &
4558 vecerr[vecerr_dim1 + 1], &errmax[nd], &nu, &nv);
4565 /* --------------------- Calculate the average error -------------
4570 mma2moy_(ndjacu, ndjacv, &ndses, &nu1, ndjacu, &nv1, ndjacv,
4571 iordru, iordrv, &patjac[idim * patjac_dim2 * patjac_dim1],
4574 /* --------------------- Set to 0.D0 the rejected coeffs ----------
4577 i__2 = idim + ndses - 1;
4578 for (ii = idim; ii <= i__2; ++ii) {
4580 for (jj = nv1; jj <= i__3; ++jj) {
4582 for (kk = nu1; kk <= i__4; ++kk) {
4583 patjac[kk + (jj + ii * patjac_dim2) * patjac_dim1] =
4592 /* --------------- Return the nb of coeff of approximation ---
4595 *ndegpu = advapp_max(*ndegpu,nu);
4596 *ndegpv = advapp_max(*ndegpv,nv);
4604 /* ------------------------------ The end -------------------------------
4606 /* --> Error in inputs */
4611 /* --------- Management of cuts, it is required 0 < NUMDEC <= 5 -------
4614 /* --> Here it is possible and necessary to cut, choose by U if it is possible */
4616 if (*numdec <= 0 || *numdec > 5) {
4625 /* --> Here it is possible and necessary to cut, choose by U if it is possible */
4627 if (*numdec <= 0 || *numdec > 5) {
4636 /* --> Here it is possible and necessary to cut, choose by 4 if it is possible */
4638 if (*numdec <= 0 || *numdec > 5) {
4643 } else if (*numdec == 2 || *numdec == 4) {
4645 } else if (*numdec == 1 || *numdec == 3) {
4653 AdvApp2Var_SysBase::maermsg_("MMA2CE2", iercod, 7L);
4655 AdvApp2Var_SysBase::mgsomsg_("MMA2CE2", 7L);
4660 //=======================================================================
4661 //function : mma2cfu_
4663 //=======================================================================
4664 int mma2cfu_(integer *ndujac,
4676 /* System generated locals */
4677 integer sosotb_dim1, disotb_dim1, disotb_offset, soditb_dim1,
4678 soditb_offset, diditb_dim1, i__1, i__2;
4680 /* Local variables */
4681 static logical ldbg;
4682 static integer nptu2, nptv2, ii, jj;
4683 static doublereal bid0, bid1, bid2;
4685 /* **********************************************************************
4690 /* Calculate the terms connected to degree NDUJAC by U of the polynomial approximation */
4691 /* of function F(u,v), starting from its discretisation
4692 /* on the roots of Legendre polynom of degree */
4693 /* NBPNTU by U and NBPNTV by V. */
4697 /* FONCTION,APPROXIMATION,COEFFICIENT,POLYNOME */
4699 /* INPUT ARGUMENTSE : */
4700 /* ------------------ */
4701 /* NDUJAC: Fixed degree by U for which the terms */
4702 /* allowing to obtain the Legendre or Jacobi coeff*/
4703 /* of even or uneven degree by V are calculated. */
4704 /* NBPNTU: Degree of Legendre polynom on the roots which of */
4705 /* the coefficients of integration by U are calculated */
4706 /* by Gauss method. It is required that NBPNTU = 30, 40, 50 or 61. */
4707 /* NBPNTV: Degree of Legendre polynom on the roots which of */
4708 /* the coefficients of integration by V are calculated */
4709 /* by Gauss method. It is required that NBPNTV = 30, 40, 50 or 61. */
4710 /* SOSOTB: Table of F(ui,vj) + F(ui,-vj) + F(-ui,vj) + F(-ui,-vj) */
4711 /* with ui and vj positive roots of Legendre polynom */
4712 /* of degree NBPNTU and NBPNTV respectively. Moreover, */
4713 /* table SOSOTB(0,j) contains F(0,vj) + F(0,-vj), */
4714 /* table SOSOTB(i,0) contains F(ui,0) + F(-ui,0) and */
4715 /* SOSOTB(0,0) contains F(0,0). */
4716 /* DISOTB: Table of F(ui,vj) + F(ui,-vj) - F(-ui,vj) - F(-ui,-vj) */
4717 /* with ui and vj positive roots of Legendre polynom */
4718 /* of degree NBPNTU and NBPNTV respectively. */
4719 /* SODITB: Table of F(ui,vj) - F(ui,-vj) + F(-ui,vj) - F(-ui,-vj) */
4720 /* with ui and vj positive roots of Legendre polynom */
4721 /* of degree NBPNTU and NBPNTV respectively. */
4722 /* DIDITB: Table of F(ui,vj) - F(ui,-vj) - F(-ui,vj) + F(-ui,-vj) */
4723 /* avec ui and vj positive roots of Legendre polynom */
4724 /* of degree NBPNTU and NBPNTV respectively. Moreover, */
4725 /* table DIDITB(0,j) contains F(0,vj) - F(0,-vj), */
4726 /* and table DIDITB(i,0) contains F(ui,0) - F(-ui,0). */
4727 /* GSSUTB: Table of coefficients of integration by Gauss method */
4728 /* Gauss by U for fixed NDUJAC : i varies from 0 to NBPNTU/2. */
4730 /* OUTPUT ARGUMENTS : */
4731 /* ------------------- */
4732 /* CHPAIR: Table of terms connected to degree NDUJAC by U to calculate the */
4733 /* coeff. of the approximation of EVEN degree by V. */
4734 /* CHIMPR: Table of terms connected to degree NDUJAC by U to calculate */
4735 /* the coeff. of approximation of UNEVEN degree by V. */
4737 /* COMMONS USED : */
4738 /* ---------------- */
4740 /* REFERENCES CALLED : */
4741 /* ----------------------- */
4743 /* DESCRIPTION/NOTES/LIMITATIONS : */
4744 /* ----------------------------------- */
4748 /* **********************************************************************
4750 /* Name of the routine */
4753 /* --------------------------- Initialisations --------------------------
4756 /* Parameter adjustments */
4758 diditb_dim1 = *nbpntu / 2 + 1;
4759 soditb_dim1 = *nbpntu / 2;
4760 soditb_offset = soditb_dim1 + 1;
4761 soditb -= soditb_offset;
4762 disotb_dim1 = *nbpntu / 2;
4763 disotb_offset = disotb_dim1 + 1;
4764 disotb -= disotb_offset;
4765 sosotb_dim1 = *nbpntu / 2 + 1;
4768 ldbg = AdvApp2Var_SysBase::mnfndeb_() >= 3;
4770 AdvApp2Var_SysBase::mgenmsg_("MMA2CFU", 7L);
4773 nptu2 = *nbpntu / 2;
4774 nptv2 = *nbpntv / 2;
4776 /* **********************************************************************
4778 /* CALCULATE COEFFICIENTS BY U */
4780 /* ----------------- Calculate coefficients of even degree --------------
4783 if (*ndujac % 2 == 0) {
4785 for (jj = 1; jj <= i__1; ++jj) {
4789 for (ii = 1; ii <= i__2; ++ii) {
4791 bid1 += sosotb[ii + jj * sosotb_dim1] * bid0;
4792 bid2 += soditb[ii + jj * soditb_dim1] * bid0;
4800 /* --------------- Calculate coefficients of uneven degree ----------
4805 for (jj = 1; jj <= i__1; ++jj) {
4809 for (ii = 1; ii <= i__2; ++ii) {
4811 bid1 += disotb[ii + jj * disotb_dim1] * bid0;
4812 bid2 += diditb[ii + jj * diditb_dim1] * bid0;
4821 /* ------- Add terms connected to the supplementary root (0.D0) ------
4822 /* ----------- of Legendre polynom of uneven degree NBPNTU -----------
4824 /* --> Only even NDUJAC terms are modified as GSSUTB(0) = 0 */
4825 /* when NDUJAC is uneven. */
4827 if (*nbpntu % 2 != 0 && *ndujac % 2 == 0) {
4830 for (jj = 1; jj <= i__1; ++jj) {
4831 chpair[jj] += sosotb[jj * sosotb_dim1] * bid0;
4832 chimpr[jj] += diditb[jj * diditb_dim1] * bid0;
4837 /* ------ Calculate the terms connected to supplementary roots (0.D0) ------
4839 /* ----------- of Legendre polynom of uneven degree NBPNTV -----------
4842 if (*nbpntv % 2 != 0) {
4843 /* --> Only CHPAIR terms are calculated as GSSVTB(0,IH-IDEBV)=0
4845 /* when IH is uneven (see MMA2CFV). */
4847 if (*ndujac % 2 == 0) {
4850 for (ii = 1; ii <= i__1; ++ii) {
4851 bid1 += sosotb[ii] * gssutb[ii];
4858 for (ii = 1; ii <= i__1; ++ii) {
4859 bid1 += diditb[ii] * gssutb[ii];
4864 if (*nbpntu % 2 != 0) {
4865 chpair[0] += sosotb[0] * gssutb[0];
4869 /* ------------------------------ The end -------------------------------
4873 AdvApp2Var_SysBase::mgsomsg_("MMA2CFU", 7L);
4878 //=======================================================================
4879 //function : mma2cfv_
4881 //=======================================================================
4882 int mma2cfv_(integer *ndvjac,
4892 /* System generated locals */
4893 integer chpair_dim1, chpair_offset, chimpr_dim1, chimpr_offset,
4894 patjac_offset, i__1, i__2;
4896 /* Local variables */
4897 static logical ldbg;
4898 static integer nptv2, ii, jj;
4899 static doublereal bid1;
4901 /* **********************************************************************
4906 /* Calculate the coefficients of polynomial approximation of F(u,v)
4907 /* of degree NDVJAC by V and of degree by U varying from MINDGU to MAXDGU.
4912 /* FONCTION,APPROXIMATION,COEFFICIENT,POLYNOME */
4914 /* INPUT ARGUMENTS : */
4915 /* ------------------ */
4917 /* NDVJAC: Degree of the polynom of approximation by V. */
4918 /* The representation in the orthogonal base starts from degre 0.
4919 /* The polynomial base is the base of Jacobi of order -1 */
4920 /* (Legendre), 0, 1 or 2 */
4921 /* MINDGU: Degree minimum by U of coeff. to calculate. */
4922 /* MAXDGU: Degree maximum by U of coeff. to calculate. */
4923 /* NBPNTV: Degree of the Legendre polynom on the roots which of */
4924 /* the coefficients of integration by V are calculated */
4925 /* by Gauss method. It is reqired that NBPNTV = 30, 40, 50 or 61 and NDVJAC < NBPNTV. */
4926 /* GSSVTB: Table of coefficients of integration by Gauss method */
4927 /* by V for NDVJAC fixed: j varies from 0 to NBPNTV/2. */
4928 /* CHPAIR: Table of terms connected to degrees from MINDGU to MAXDGU by U to
4929 /* calculate the coeff. of approximation of EVEN degree NDVJAC by V. */
4930 /* CHIMPR: Table of terms connected to degrees from MINDGU to MAXDGU by U to
4931 /* calculate the coeff. of approximation of UNEVEN degree NDVJAC by V. */
4933 /* OUTPUT ARGUMENTS : */
4934 /* ------------------- */
4935 /* PATJAC: Table of coefficients by U of the polynom of approximation */
4936 /* P(u,v) of degree MINDGU to MAXDGU by U and NDVJAC by V. */
4938 /* COMMONS USED : */
4939 /* -------------- */
4941 /* REFERENCES CALLED : */
4942 /* --------------------- */
4944 /* DESCRIPTION/NOTES/LIMITATIONS : */
4945 /* ------------------------------- */
4947 /* **********************************************************************
4949 /* Name of the routine */
4952 /* --------------------------- Initialisations --------------------------
4955 /* Parameter adjustments */
4956 patjac_offset = *mindgu;
4957 patjac -= patjac_offset;
4958 chimpr_dim1 = *nbpntv / 2;
4959 chimpr_offset = chimpr_dim1 * *mindgu + 1;
4960 chimpr -= chimpr_offset;
4961 chpair_dim1 = *nbpntv / 2 + 1;
4962 chpair_offset = chpair_dim1 * *mindgu;
4963 chpair -= chpair_offset;
4966 ldbg = AdvApp2Var_SysBase::mnfndeb_() >= 3;
4968 AdvApp2Var_SysBase::mgenmsg_("MMA2CFV", 7L);
4970 nptv2 = *nbpntv / 2;
4972 /* --------- Calculate the coefficients for even degree NDVJAC ----------
4975 if (*ndvjac % 2 == 0) {
4977 for (ii = *mindgu; ii <= i__1; ++ii) {
4980 for (jj = 1; jj <= i__2; ++jj) {
4981 bid1 += chpair[jj + ii * chpair_dim1] * gssvtb[jj];
4988 /* -------- Calculate the coefficients for uneven degree NDVJAC -----
4993 for (ii = *mindgu; ii <= i__1; ++ii) {
4996 for (jj = 1; jj <= i__2; ++jj) {
4997 bid1 += chimpr[jj + ii * chimpr_dim1] * gssvtb[jj];
5005 /* ------- Add terms connected to the supplementary root (0.D0) ----- */
5006 /* --------of the Legendre polynom of uneven degree NBPNTV --------- */
5008 if (*nbpntv % 2 != 0 && *ndvjac % 2 == 0) {
5011 for (ii = *mindgu; ii <= i__1; ++ii) {
5012 patjac[ii] += bid1 * chpair[ii * chpair_dim1];
5017 /* ------------------------------ The end -------------------------------
5021 AdvApp2Var_SysBase::mgsomsg_("MMA2CFV", 7L);
5026 //=======================================================================
5027 //function : mma2ds1_
5029 //=======================================================================
5030 int AdvApp2Var_ApproxF2var::mma2ds1_(integer *ndimen,
5033 const AdvApp2Var_EvaluatorFunc2Var& foncnp,
5048 /* System generated locals */
5049 integer sosotb_dim1, sosotb_dim2, sosotb_offset, disotb_dim1, disotb_dim2,
5050 disotb_offset, soditb_dim1, soditb_dim2, soditb_offset,
5051 diditb_dim1, diditb_dim2, diditb_offset, fpntab_dim1,
5052 fpntab_offset, i__1;
5054 /* Local variables */
5055 static logical ldbg;
5056 static integer ibid1, ibid2, iuouv, nd;
5057 static integer isz1, isz2;
5059 /* **********************************************************************
5064 /* Discretisation of function F(u,v) on the roots of Legendre polynoms. */
5068 /* FONCTION&,DISCRETISATION,&POINT */
5070 /* INPUT ARGUMENTS : */
5071 /* ------------------ */
5072 /* NDIMEN: Dimension of the space. */
5073 /* UINTFN: Limits of the interval of definition by u of the function */
5074 /* to be processed: (UINTFN(1),UINTFN(2)). */
5075 /* VINTFN: Limits of the interval of definition by v of the function */
5076 /* to be processed: (VINTFN(1),VINTFN(2)). */
5077 /* FONCNP: The NAME of the non-polynomial function to be processed. */
5078 /* NBPNTU: The degree of Legendre polynom on the roots which of */
5079 /* FONCNP is discretized by u. */
5080 /* NBPNTV: The degree of Legendre polynom on the roots which of */
5081 /* FONCNP is discretized by v. */
5082 /* UROOTB: Table of STRICTLY POSITIVE roots of the polynom */
5083 /* of Legendre of degree NBPNTU defined on (-1,1). */
5084 /* VROOTB: Table of STRICTLY POSITIVE roots of the polynom */
5085 /* of Legendre of degree NBPNTV defined on (-1,1). */
5086 /* ISOFAV: Shows the type of iso of F(u,v) to be extracted to improve */
5087 /* the rapidity of calculation (has no influence on the form */
5089 /* = 1, shows that it is necessary to calculate the points of F(u,v) */
5090 /* with fixed u (with NBPNTV values different from v). */
5091 /* = 2, shows that it is necessaty to calculate the points of F(u,v) */
5092 /* with fixed v (with NBPNTU values different from u). */
5093 /* SOSOTB: Preinitialized table (input/output argument). */
5094 /* DISOTB: Preinitialized table (input/output argument). */
5095 /* SODITB: Preinitialized table (input/output argument). */
5096 /* DIDITB: Preinitialized table (input/output argument). */
5098 /* OUTPUT ARGUMENTS : */
5099 /* ------------------- */
5100 /* SOSOTB: Table where the terms */
5101 /* F(ui,vj) + F(ui,-vj) + F(-ui,vj) + F(-ui,-vj) */
5102 /* are added with ui and vj positive roots of Legendre polynom */
5103 /* of degree NBPNTU and NBPNTV respectively. */
5104 /* DISOTB: Table where the terms */
5105 /* F(ui,vj) + F(ui,-vj) - F(-ui,vj) - F(-ui,-vj) */
5106 /* are added with ui and vj positive roots of Legendre polynom */
5107 /* of degree NBPNTU and NBPNTV respectively. */
5108 /* SODITB: Table where the terms */
5109 /* F(ui,vj) - F(ui,-vj) + F(-ui,vj) - F(-ui,-vj) */
5110 /* are added with ui and vj positive roots of Legendre polynom */
5111 /* of degree NBPNTU and NBPNTV respectively. */
5112 /* DIDITB: Table where the terms */
5113 /* F(ui,vj) - F(ui,-vj) - F(-ui,vj) + F(-ui,-vj) */
5114 /* are added with ui and vj positive roots of Legendre polynom */
5115 /* of degree NBPNTU and NBPNTV respectively. */
5116 /* FPNTAB: Auxiliary table. */
5117 /* TTABLE: Auxiliary table. */
5118 /* IERCOD: Error code >100 Pb in the evaluation of FONCNP, */
5119 /* the returned error code is equal to error code of FONCNP + 100. */
5121 /* COMMONS USED : */
5122 /* ---------------- */
5124 /* REFERENCES CALLED : */
5125 /* --------------------- */
5127 /* DESCRIPTION/NOTES/LIMITATIONS : */
5128 /* ----------------------------------- */
5129 /* --> The external function created by the caller of MA2F1K, MA2FDK */
5130 /* where MA2FXK should be in the following form : */
5131 /* SUBROUTINE FONCNP(NDIMEN,UINTFN,VINTFN,ISOFAV,TCONST,NBPTAB */
5132 /* ,TTABLE,IDERIU,IDERIV,PPNTAB,IERCOD) */
5133 /* with the following input arguments : */
5134 /* - NDIMEN is integer defined as the sum of dimensions of */
5135 /* sub-spaces (i.e. total dimension of the problem). */
5136 /* - UINTFN(2) is a table of 2 reals containing the interval */
5137 /* by u where the function to be approximated is defined */
5138 /* (so it is equal to UIFONC). */
5139 /* - VINTFN(2) is a table of 2 reals containing the interval */
5140 /* by v where the function to be approximated is defined */
5141 /* (so it is equal to VIFONC). */
5142 /* - ISOFAV, is 1 if it is necessary to calculate points with constant u, */
5143 /* is 2 if it is necessary to calculate points with constant v. */
5144 /* Any other value is an error. */
5145 /* - TCONST, real, value of the fixed parameter. Takes values */
5146 /* in (UIFONC(1),UIFONC(2)) if ISOFAV = 1 or */
5147 /* ins (VIFONC(1),VIFONC(2)) if ISOFAV = 2. */
5148 /* - NBPTAB, integer. Shows the number of points to be calculated. */
5149 /* - TTABLE, a table of reals NBPTAB. These are the values of */
5150 /* 'free' parameter of discretization (v if IISOFAV=1, */
5151 /* u if IISOFAV=2). */
5152 /* - IDERIU, integer, takes values between 0 (position) */
5153 /* and IORDRE(1) (partial derivative of the function by u */
5154 /* of order IORDRE(1) if IORDRE(1) > 0). */
5155 /* - IDERIV, integer, takes values between 0 (position) */
5156 /* and IORDRE(2) (partial derivative of the function by v */
5157 /* of order IORDRE(2) if IORDRE(2) > 0). */
5158 /* If IDERIU=i and IDERIV=j, FONCNP should calculate the */
5159 /* points of the derivative : */
5166 /* and the output arguments aret : */
5167 /* - FPNTAB(NDIMEN,NBPTAB) contains, at output, the table of */
5168 /* NBPTAB points calculated in FONCNP. */
5169 /* - IERCOD is, at output the error code of FONCNP. This code */
5170 /* (integer) should be strictly positive if there is a problem. */
5172 /* The input arguments SHOULD NOT be modified under FONCNP.
5175 /* -->As FONCNP is not forcedly defined in (-1,1)*(-1,1), the */
5176 /* values of UROOTB and VROOTB are consequently modified. */
5178 /* -->The results of discretisation are ranked in 4 tables */
5179 /* SOSOTB, DISOTB, SODITB and DIDITB to earn time */
5180 /* during the calculation of coefficients of the polynom of approximation. */
5182 /* When NBPNTU is uneven : */
5183 /* table SOSOTB(0,j) contains F(0,vj) + F(0,-vj), */
5184 /* table DIDITB(0,j) contains F(0,vj) - F(0,-vj), */
5185 /* When NBPNTV is uneven : */
5186 /* table SOSOTB(i,0) contains F(ui,0) + F(-ui,0), */
5187 /* table DIDITB(i,0) contains F(ui,0) - F(-ui,0), */
5188 /* When NBPNTU and NBPNTV are uneven : */
5189 /* term SOSOTB(0,0) contains F(0,0). */
5192 /* **********************************************************************
5194 /* Name of the routine */
5197 /* --------------------------- Initialization --------------------------
5200 /* Parameter adjustments */
5201 fpntab_dim1 = *ndimen;
5202 fpntab_offset = fpntab_dim1 + 1;
5203 fpntab -= fpntab_offset;
5207 diditb_dim1 = *nbpntu / 2 + 1;
5208 diditb_dim2 = *nbpntv / 2 + 1;
5209 diditb_offset = diditb_dim1 * diditb_dim2;
5210 diditb -= diditb_offset;
5211 soditb_dim1 = *nbpntu / 2;
5212 soditb_dim2 = *nbpntv / 2;
5213 soditb_offset = soditb_dim1 * (soditb_dim2 + 1) + 1;
5214 soditb -= soditb_offset;
5215 disotb_dim1 = *nbpntu / 2;
5216 disotb_dim2 = *nbpntv / 2;
5217 disotb_offset = disotb_dim1 * (disotb_dim2 + 1) + 1;
5218 disotb -= disotb_offset;
5219 sosotb_dim1 = *nbpntu / 2 + 1;
5220 sosotb_dim2 = *nbpntv / 2 + 1;
5221 sosotb_offset = sosotb_dim1 * sosotb_dim2;
5222 sosotb -= sosotb_offset;
5227 ldbg = AdvApp2Var_SysBase::mnfndeb_() >= 3;
5229 AdvApp2Var_SysBase::mgenmsg_("MMA2DS1", 7L);
5232 if (*isofav < 1 || *isofav > 2) {
5238 /* **********************************************************************
5240 /* --------- Discretization by U on the roots of the polynom of ------ */
5241 /* --------------- Legendre of degree NBPNTU, iso-V by iso-V --------- */
5242 /* **********************************************************************
5246 mma2ds2_(ndimen, &uintfn[1], &vintfn[1], foncnp, nbpntu, nbpntv, &
5247 urootb[1], &vrootb[1], &iuouv, &sosotb[sosotb_offset], &
5248 disotb[disotb_offset], &soditb[soditb_offset], &diditb[
5249 diditb_offset], &fpntab[fpntab_offset], &ttable[1], iercod);
5251 /* ******************************************************************
5253 /* --------- Discretization by V on the roots of the polynom of ------ */
5254 /* --------------- Legendre of degree NBPNTV, iso-V by iso-V --------- */
5255 /* ******************************************************************
5259 /* --> Inversion of indices of tables */
5261 for (nd = 1; nd <= i__1; ++nd) {
5262 isz1 = *nbpntu / 2 + 1;
5263 isz2 = *nbpntv / 2 + 1;
5264 AdvApp2Var_MathBase::mmfmtb1_(&isz1, &sosotb[nd * sosotb_dim2 * sosotb_dim1], &isz1, &
5265 isz2, &isz2, &sosotb[nd * sosotb_dim2 * sosotb_dim1], &
5266 ibid1, &ibid2, iercod);
5270 AdvApp2Var_MathBase::mmfmtb1_(&isz1, &diditb[nd * diditb_dim2 * diditb_dim1], &isz1, &
5271 isz2, &isz2, &diditb[nd * diditb_dim2 * diditb_dim1], &
5272 ibid1, &ibid2, iercod);
5278 AdvApp2Var_MathBase::mmfmtb1_(&isz1, &soditb[(nd * soditb_dim2 + 1) * soditb_dim1 + 1],
5279 &isz1, &isz2, &isz2, &soditb[(nd * soditb_dim2 + 1) *
5280 soditb_dim1 + 1], &ibid1, &ibid2, iercod);
5284 AdvApp2Var_MathBase::mmfmtb1_(&isz1, &disotb[(nd * disotb_dim2 + 1) * disotb_dim1 + 1],
5285 &isz1, &isz2, &isz2, &disotb[(nd * disotb_dim2 + 1) *
5286 disotb_dim1 + 1], &ibid1, &ibid2, iercod);
5293 mma2ds2_(ndimen, &vintfn[1], &uintfn[1], foncnp, nbpntv, nbpntu, &
5294 vrootb[1], &urootb[1], &iuouv, &sosotb[sosotb_offset], &
5295 soditb[soditb_offset], &disotb[disotb_offset], &diditb[
5296 diditb_offset], &fpntab[fpntab_offset], &ttable[1], iercod);
5297 /* --> Inversion of indices of tables */
5299 for (nd = 1; nd <= i__1; ++nd) {
5300 isz1 = *nbpntv / 2 + 1;
5301 isz2 = *nbpntu / 2 + 1;
5302 AdvApp2Var_MathBase::mmfmtb1_(&isz1, &sosotb[nd * sosotb_dim2 * sosotb_dim1], &isz1, &
5303 isz2, &isz2, &sosotb[nd * sosotb_dim2 * sosotb_dim1], &
5304 ibid1, &ibid2, iercod);
5308 AdvApp2Var_MathBase::mmfmtb1_(&isz1, &diditb[nd * diditb_dim2 * diditb_dim1], &isz1, &
5309 isz2, &isz2, &diditb[nd * diditb_dim2 * diditb_dim1], &
5310 ibid1, &ibid2, iercod);
5316 AdvApp2Var_MathBase::mmfmtb1_(&isz1, &soditb[(nd * soditb_dim2 + 1) * soditb_dim1 + 1],
5317 &isz1, &isz2, &isz2, &soditb[(nd * soditb_dim2 + 1) *
5318 soditb_dim1 + 1], &ibid1, &ibid2, iercod);
5322 AdvApp2Var_MathBase::mmfmtb1_(&isz1, &disotb[(nd * disotb_dim2 + 1) * disotb_dim1 + 1],
5323 &isz1, &isz2, &isz2, &disotb[(nd * disotb_dim2 + 1) *
5324 disotb_dim1 + 1], &ibid1, &ibid2, iercod);
5332 /* ------------------------------ The end -------------------------------
5338 AdvApp2Var_SysBase::maermsg_("MMA2DS1", iercod, 7L);
5341 AdvApp2Var_SysBase::mgsomsg_("MMA2DS1", 7L);
5346 //=======================================================================
5347 //function : mma2ds2_
5349 //=======================================================================
5350 int mma2ds2_(integer *ndimen,
5353 const AdvApp2Var_EvaluatorFunc2Var& foncnp,
5368 static integer c__0 = 0;
5369 /* System generated locals */
5370 integer sosotb_dim1, sosotb_dim2, sosotb_offset, disotb_dim1, disotb_dim2,
5371 disotb_offset, soditb_dim1, soditb_dim2, soditb_offset,
5372 diditb_dim1, diditb_dim2, diditb_offset, fpntab_dim1,
5373 fpntab_offset, i__1, i__2, i__3;
5375 /* Local variables */
5376 static integer jdec;
5377 static logical ldbg;
5378 static doublereal alinu, blinu, alinv, blinv, tcons;
5379 static doublereal dbfn1[2], dbfn2[2];
5380 static integer nuroo, nvroo, id, iu, iv;
5381 static doublereal um, up;
5384 /* **********************************************************************
5389 /* Discretization of function F(u,v) on the roots of polynoms of Legendre. */
5393 /* FONCTION&,DISCRETISATION,&POINT */
5395 /* INPUT ARGUMENTS : */
5396 /* ------------------ */
5397 /* NDIMEN: Dimension of the space. */
5398 /* UINTFN: Limits of the interval of definition by u of the function */
5399 /* to be processed: (UINTFN(1),UINTFN(2)). */
5400 /* VINTFN: Limits of the interval of definition by v of the function */
5401 /* to be processed: (VINTFN(1),VINTFN(2)). */
5402 /* FONCNP: The NAME of the non-polynomial function to be processed. */
5403 /* NBPNTU: The degree of Legendre polynom on the roots which of */
5404 /* FONCNP is discretized by u. */
5405 /* NBPNTV: The degree of Legendre polynom on the roots which of */
5406 /* FONCNP is discretized by v. */
5407 /* UROOTB: Table of STRICTLY POSITIVE roots of the polynom */
5408 /* of Legendre of degree NBPNTU defined on (-1,1). */
5409 /* VROOTB: Table of STRICTLY POSITIVE roots of the polynom */
5410 /* of Legendre of degree NBPNTV defined on (-1,1). */
5411 /* IIUOUV: Shows the type of iso of F(u,v) tom be extracted to improve the */
5412 /* rapidity of calculation (has no influence on the form of result) */
5413 /* = 1, shows that it is necessary to calculate the points of F(u,v) */
5414 /* with fixed u (so with NBPNTV values different from v). */
5415 /* = 2, shows that it is necessary to calculate the points of F(u,v) */
5416 /* with fixed v (so with NBPNTV values different from u). */
5417 /* SOSOTB: Preinitialized table (input/output argument). */
5418 /* DISOTB: Preinitialized table (input/output argument). */
5419 /* SODITB: Preinitialized table (input/output argument). */
5420 /* DIDITB: Preinitialized table (input/output argument). */
5422 /* OUTPUT ARGUMENTS : */
5423 /* ------------------- */
5424 /* SOSOTB: Table where the terms */
5425 /* F(ui,vj) + F(ui,-vj) + F(-ui,vj) + F(-ui,-vj) */
5426 /* are added with ui and vj positive roots of Legendre polynom */
5427 /* of degree NBPNTU and NBPNTV respectively. */
5428 /* DISOTB: Table where the terms */
5429 /* F(ui,vj) + F(ui,-vj) - F(-ui,vj) - F(-ui,-vj) */
5430 /* are added with ui and vj positive roots of Legendre polynom */
5431 /* of degree NBPNTU and NBPNTV respectively. */
5432 /* SODITB: Table where the terms */
5433 /* F(ui,vj) - F(ui,-vj) + F(-ui,vj) - F(-ui,-vj) */
5434 /* are added with ui and vj positive roots of Legendre polynom */
5435 /* of degree NBPNTU and NBPNTV respectively. */
5436 /* DIDITB: Table where the terms */
5437 /* F(ui,vj) - F(ui,-vj) - F(-ui,vj) + F(-ui,-vj) */
5438 /* are added with ui and vj positive roots of Legendre polynom */
5439 /* of degree NBPNTU and NBPNTV respectively. */
5440 /* FPNTAB: Auxiliary table. */
5441 /* TTABLE: Auxiliary table. */
5442 /* IERCOD: Error code >100 Pb in the evaluation of FONCNP, */
5443 /* the returned error code is equal to error code of FONCNP + 100. */
5445 /* COMMONS USED : */
5446 /* ---------------- */
5448 /* REFERENCES CALLED : */
5449 /* --------------------- */
5451 /* DESCRIPTION/NOTES/LIMITATIONS : */
5452 /* ----------------------------------- */
5453 /* --> The external function created by the caller of MA2F1K, MA2FDK */
5454 /* where MA2FXK should be in the following form : */
5455 /* SUBROUTINE FONCNP(NDIMEN,UINTFN,VINTFN,IIIUOUV,TCONST,NBPTAB */
5456 /* ,TTABLE,IDERIU,IDERIV,PPNTAB,IERCOD) */
5457 /* with the following input arguments : */
5458 /* - NDIMEN is integer defined as the sum of dimensions of */
5459 /* sub-spaces (i.e. total dimension of the problem). */
5460 /* - UINTFN(2) is a table of 2 reals containing the interval */
5461 /* by u where the function to be approximated is defined */
5462 /* (so it is equal to UIFONC). */
5463 /* - VINTFN(2) is a table of 2 reals containing the interval */
5464 /* by v where the function to be approximated is defined */
5465 /* (so it is equal to VIFONC). */
5466 /* - IIIUOUV, is 1 if it is necessary to calculate points with constant u, */
5467 /* is 2 if it is necessary to calculate points with constant v. */
5468 /* Any other value is an error. */
5469 /* - TCONST, real, value of the fixed parameter. Takes values */
5470 /* in (UIFONC(1),UIFONC(2)) if ISOFAV = 1 or */
5471 /* ins (VIFONC(1),VIFONC(2)) if ISOFAV = 2. */
5472 /* - NBPTAB, integer. Shows the number of points to be calculated. */
5473 /* - TTABLE, a table of reals NBPTAB. These are the values of */
5474 /* 'free' parameter of discretization (v if IIIUOUV=1, */
5475 /* u if IIIUOUV=2). */
5476 /* - IDERIU, integer, takes values between 0 (position) */
5477 /* and IORDRE(1) (partial derivative of the function by u */
5478 /* of order IORDRE(1) if IORDRE(1) > 0). */
5479 /* - IDERIV, integer, takes values between 0 (position) */
5480 /* and IORDRE(2) (partial derivative of the function by v */
5481 /* of order IORDRE(2) if IORDRE(2) > 0). */
5482 /* If IDERIU=i and IDERIV=j, FONCNP should calculate the */
5483 /* points of the derivative : */
5490 /* and the output arguments aret : */
5491 /* - FPNTAB(NDIMEN,NBPTAB) contains, at output, the table of */
5492 /* NBPTAB points calculated in FONCNP. */
5493 /* - IERCOD is, at output the error code of FONCNP. This code */
5494 /* (integer) should be strictly positive if there is a problem. */
5496 /* The input arguments SHOULD NOT be modified under FONCNP.
5499 /* -->As FONCNP is not forcedly defined in (-1,1)*(-1,1), the */
5500 /* values of UROOTB and VROOTB are consequently modified. */
5502 /* -->The results of discretisation are ranked in 4 tables */
5503 /* SOSOTB, DISOTB, SODITB and DIDITB to earn time */
5504 /* during the calculation of coefficients of the polynom of approximation. */
5506 /* When NBPNTU is uneven : */
5507 /* table SOSOTB(0,j) contains F(0,vj) + F(0,-vj), */
5508 /* table DIDITB(0,j) contains F(0,vj) - F(0,-vj), */
5509 /* When NBPNTV is uneven : */
5510 /* table SOSOTB(i,0) contains F(ui,0) + F(-ui,0), */
5511 /* table DIDITB(i,0) contains F(ui,0) - F(-ui,0), */
5512 /* When NBPNTU and NBPNTV are uneven : */
5513 /* term SOSOTB(0,0) contains F(0,0). */
5515 /* ATTENTION: These 4 tables are filled by varying the */
5516 /* 1st index first. So, the discretizations */
5517 /* of F(...,t) (for IIUOUV = 2) or of F(t,...) (IIUOUV = 1) */
5518 /* are stored in SOSOTB(...,t), SODITB(...,t), etc... */
5519 /* (this allows to gain important time). */
5520 /* It is required that the caller, in case of IIUOUV=1, */
5521 /* invert the roles of u and v, of SODITB and DISOTB BEFORE the */
5524 /* **********************************************************************
5527 /* Name of the routine */
5529 /* --> Indices of loops. */
5531 /* --------------------------- Initialization --------------------------
5534 /* Parameter adjustments */
5538 fpntab_dim1 = *ndimen;
5539 fpntab_offset = fpntab_dim1 + 1;
5540 fpntab -= fpntab_offset;
5542 diditb_dim1 = *nbpntu / 2 + 1;
5543 diditb_dim2 = *nbpntv / 2 + 1;
5544 diditb_offset = diditb_dim1 * diditb_dim2;
5545 diditb -= diditb_offset;
5546 soditb_dim1 = *nbpntu / 2;
5547 soditb_dim2 = *nbpntv / 2;
5548 soditb_offset = soditb_dim1 * (soditb_dim2 + 1) + 1;
5549 soditb -= soditb_offset;
5550 disotb_dim1 = *nbpntu / 2;
5551 disotb_dim2 = *nbpntv / 2;
5552 disotb_offset = disotb_dim1 * (disotb_dim2 + 1) + 1;
5553 disotb -= disotb_offset;
5554 sosotb_dim1 = *nbpntu / 2 + 1;
5555 sosotb_dim2 = *nbpntv / 2 + 1;
5556 sosotb_offset = sosotb_dim1 * sosotb_dim2;
5557 sosotb -= sosotb_offset;
5561 ldbg = AdvApp2Var_SysBase::mnfndeb_() >= 3;
5563 AdvApp2Var_SysBase::mgenmsg_("MMA2DS2", 7L);
5567 alinu = (uintfn[2] - uintfn[1]) / 2.;
5568 blinu = (uintfn[2] + uintfn[1]) / 2.;
5569 alinv = (vintfn[2] - vintfn[1]) / 2.;
5570 blinv = (vintfn[2] + vintfn[1]) / 2.;
5573 dbfn1[0] = vintfn[1];
5574 dbfn1[1] = vintfn[2];
5575 dbfn2[0] = uintfn[1];
5576 dbfn2[1] = uintfn[2];
5578 dbfn1[0] = uintfn[1];
5579 dbfn1[1] = uintfn[2];
5580 dbfn2[0] = vintfn[1];
5581 dbfn2[1] = vintfn[2];
5584 /* **********************************************************************
5586 /* -------- Discretization by U on the roots of Legendre polynom -------- */
5587 /* ---------------- of degree NBPNTU, with Vj fixed -------------------- */
5588 /* **********************************************************************
5591 nuroo = *nbpntu / 2;
5592 nvroo = *nbpntv / 2;
5593 jdec = (*nbpntu + 1) / 2;
5595 /* ----------- Loading of parameters of discretization by U ------------- */
5598 for (iu = 1; iu <= i__1; ++iu) {
5599 ttable[iu] = blinu + alinu * urootb[iu];
5603 /* -------------- For Vj fixed, negative root of Legendre ------------- */
5606 for (iv = 1; iv <= i__1; ++iv) {
5607 tcons = blinv + alinv * vrootb[iv];
5608 (*const_cast <AdvApp2Var_EvaluatorFunc2Var*> (&foncnp)).Evaluate (
5609 ndimen, dbfn1, dbfn2, iiuouv, &tcons, nbpntu, &
5610 ttable[1], &c__0, &c__0, &fpntab[fpntab_offset], iercod);
5615 for (id = 1; id <= i__2; ++id) {
5617 for (iu = 1; iu <= i__3; ++iu) {
5618 up = fpntab[id + (iu + jdec) * fpntab_dim1];
5619 um = fpntab[id + (nuroo - iu + 1) * fpntab_dim1];
5620 sosotb[iu + (nvroo - iv + 1 + id * sosotb_dim2) * sosotb_dim1]
5621 = sosotb[iu + (nvroo - iv + 1 + id * sosotb_dim2) *
5622 sosotb_dim1] + up + um;
5623 disotb[iu + (nvroo - iv + 1 + id * disotb_dim2) * disotb_dim1]
5624 = disotb[iu + (nvroo - iv + 1 + id * disotb_dim2) *
5625 disotb_dim1] + up - um;
5626 soditb[iu + (nvroo - iv + 1 + id * soditb_dim2) * soditb_dim1]
5627 = soditb[iu + (nvroo - iv + 1 + id * soditb_dim2) *
5628 soditb_dim1] - up - um;
5629 diditb[iu + (nvroo - iv + 1 + id * diditb_dim2) * diditb_dim1]
5630 = diditb[iu + (nvroo - iv + 1 + id * diditb_dim2) *
5631 diditb_dim1] - up + um;
5634 if (*nbpntu % 2 != 0) {
5635 up = fpntab[id + jdec * fpntab_dim1];
5636 sosotb[(nvroo - iv + 1 + id * sosotb_dim2) * sosotb_dim1] +=
5638 diditb[(nvroo - iv + 1 + id * diditb_dim2) * diditb_dim1] -=
5646 /* --------- For Vj = 0 (uneven NBPNTV), discretization by U ----------- */
5648 if (*nbpntv % 2 != 0) {
5650 (*const_cast <AdvApp2Var_EvaluatorFunc2Var*> (&foncnp)).Evaluate (
5651 ndimen, dbfn1, dbfn2, iiuouv, &tcons, nbpntu, &
5652 ttable[1], &c__0, &c__0, &fpntab[fpntab_offset], iercod);
5657 for (id = 1; id <= i__1; ++id) {
5659 for (iu = 1; iu <= i__2; ++iu) {
5660 up = fpntab[id + (jdec + iu) * fpntab_dim1];
5661 um = fpntab[id + (nuroo - iu + 1) * fpntab_dim1];
5662 sosotb[iu + id * sosotb_dim2 * sosotb_dim1] = sosotb[iu + id *
5663 sosotb_dim2 * sosotb_dim1] + up + um;
5664 diditb[iu + id * diditb_dim2 * diditb_dim1] = diditb[iu + id *
5665 diditb_dim2 * diditb_dim1] + up - um;
5668 if (*nbpntu % 2 != 0) {
5669 up = fpntab[id + jdec * fpntab_dim1];
5670 sosotb[id * sosotb_dim2 * sosotb_dim1] += up;
5676 /* -------------- For Vj fixed, positive root of Legendre ------------- */
5679 for (iv = 1; iv <= i__1; ++iv) {
5680 tcons = alinv * vrootb[(*nbpntv + 1) / 2 + iv] + blinv;
5681 (*const_cast <AdvApp2Var_EvaluatorFunc2Var*> (&foncnp)).Evaluate (
5682 ndimen, dbfn1, dbfn2, iiuouv, &tcons, nbpntu, &
5683 ttable[1], &c__0, &c__0, &fpntab[fpntab_offset], iercod);
5688 for (id = 1; id <= i__2; ++id) {
5690 for (iu = 1; iu <= i__3; ++iu) {
5691 up = fpntab[id + (iu + jdec) * fpntab_dim1];
5692 um = fpntab[id + (nuroo - iu + 1) * fpntab_dim1];
5693 sosotb[iu + (iv + id * sosotb_dim2) * sosotb_dim1] = sosotb[
5694 iu + (iv + id * sosotb_dim2) * sosotb_dim1] + up + um;
5695 disotb[iu + (iv + id * disotb_dim2) * disotb_dim1] = disotb[
5696 iu + (iv + id * disotb_dim2) * disotb_dim1] + up - um;
5697 soditb[iu + (iv + id * soditb_dim2) * soditb_dim1] = soditb[
5698 iu + (iv + id * soditb_dim2) * soditb_dim1] + up + um;
5699 diditb[iu + (iv + id * diditb_dim2) * diditb_dim1] = diditb[
5700 iu + (iv + id * diditb_dim2) * diditb_dim1] + up - um;
5703 if (*nbpntu % 2 != 0) {
5704 up = fpntab[id + jdec * fpntab_dim1];
5705 sosotb[(iv + id * sosotb_dim2) * sosotb_dim1] += up;
5706 diditb[(iv + id * diditb_dim2) * diditb_dim1] += up;
5713 /* ------------------------------ The end -------------------------------
5719 AdvApp2Var_SysBase::maermsg_("MMA2DS2", iercod, 7L);
5722 AdvApp2Var_SysBase::mgsomsg_("MMA2DS2", 7L);
5727 //=======================================================================
5728 //function : mma2er1_
5730 //=======================================================================
5731 int mma2er1_(integer *ndjacu,
5747 /* System generated locals */
5748 integer patjac_dim1, patjac_dim2, patjac_offset, i__1, i__2, i__3;
5751 /* Local variables */
5752 static logical ldbg;
5753 static integer minu, minv;
5754 static doublereal vaux[2];
5755 static integer ii, nd, jj;
5756 static doublereal bid0, bid1;
5758 /* **********************************************************************
5763 /* Calculate max approximation error done when */
5764 /* the coefficients of PATJAC such that the degree by U varies between */
5765 /* MINDGU and MAXDGU and the degree by V varies between MINDGV and MAXDGV are removed. */
5769 /* TOUS,AB_SPECIFI:: CARREAU&,CALCUL,&ERREUR */
5771 /* INPUT ARGUMENTS : */
5772 /* ------------------ */
5773 /* NDJACU: Dimension by U of table PATJAC. */
5774 /* NDJACV: Dimension by V of table PATJAC. */
5775 /* NDIMEN: Dimension of the space. */
5776 /* MINDGU: Lower limit of index by U of coeff. of PATJAC to be taken into account. */
5777 /* MAXDGU: Upper limit of index by U of coeff. of PATJAC to be taken into account. */
5778 /* MINDGV: Lower limit of index by V of coeff. of PATJAC to be taken into account. */
5779 /* MAXDGV: Upper limit of index by V of coeff. of PATJAC to be taken into account. */
5780 /* IORDRU: Order of continuity by U provided by square PATJAC (from -1 to 2) */
5781 /* IORDRV: Order of continuity by U provided by square PATJAC (from -1 to 2) */
5782 /* XMAXJU: Maximum value of Jacobi polynoms of order IORDRU, */
5783 /* from degree 0 to MAXDGU - 2*(IORDU+1) */
5784 /* XMAXJV: Maximum value of Jacobi polynoms of order IORDRV, */
5785 /* from degree 0 to MAXDGV - 2*(IORDV+1) */
5786 /* PATJAC: Table of coeff. of square of approximation with */
5787 /* constraints of order IORDRU by U and IORDRV by V. */
5788 /* VECERR: Auxiliary vector. */
5789 /* ERREUR: MAX Error commited during removal of ALREADY CALCULATED coeff of PATJAC */
5791 /* OUTPUT ARGUMENTS : */
5792 /* ------------------- */
5793 /* ERREUR: MAX Error commited during removal of coeff of PATJAC */
5794 /* of indices from MINDGU to MAXDGU by U and from MINDGV to MAXDGV by V */
5795 /* THEN the already calculated error. */
5797 /* COMMONS USED : */
5798 /* ---------------- */
5800 /* REFERENCES CALLED : */
5801 /* --------------------- */
5803 /* DESCRIPTION/NOTES/LIMITATIONS : */
5804 /* ----------------------------------- */
5805 /* Table PATJAC is the place of storage of coeff. Cij of the square of */
5806 /* approximation of F(U,V). The indices i and j show the degree */
5807 /* by U and by V of base polynoms. These polynoms have the form: */
5809 /* ((1 - U*U)**(IORDRU+1)).J(i-2*(IORDRU+1)(U), where */
5811 /* polynom J(i-2*(IORDU+1)(U) is the Jacobi polynom of order */
5812 /* IORDRU+1 (the same by V by replacing U u V in the expression above). */
5814 /* The contribution to the error of term Cij when it is */
5815 /* removed from PATJAC is increased by: */
5817 /* DABS(Cij)*XMAXJU(i-2*(IORDRU+1))*XMAXJV(J-2*(IORDRV+1)) where */
5819 /* XMAXJU(i-2*(IORDRU+1) = ((1 - U*U)**(IORDRU+1)).J(i-2*(IORDRU+1)(U),
5821 /* XMAXJV(i-2*(IORDRV+1) = ((1 - V*V)**(IORDRV+1)).J(j-2*(IORDRV+1)(V).
5825 /* ***********************************************************************
5827 /* Name of the routine */
5830 /* ----------------------------- Initialisations ------------------------
5833 /* Parameter adjustments */
5835 patjac_dim1 = *ndjacu + 1;
5836 patjac_dim2 = *ndjacv + 1;
5837 patjac_offset = patjac_dim1 * patjac_dim2;
5838 patjac -= patjac_offset;
5841 ldbg = AdvApp2Var_SysBase::mnfndeb_() >= 3;
5843 AdvApp2Var_SysBase::mgenmsg_("MMA2ER1", 7L);
5846 minu = (*iordru + 1) << 1;
5847 minv = (*iordrv + 1) << 1;
5849 /* ------------------- Calculate the increment of the max error --------------- */
5850 /* ----- during the removal of the coeffs of indices from MINDGU to MAXDGU ---- */
5851 /* ---------------- by U and indices from MINDGV to MAXDGV by V --------------- */
5854 for (nd = 1; nd <= i__1; ++nd) {
5857 for (jj = *mindgv; jj <= i__2; ++jj) {
5860 for (ii = *mindgu; ii <= i__3; ++ii) {
5861 bid0 += (d__1 = patjac[ii + (jj + nd * patjac_dim2) *
5862 patjac_dim1], advapp_abs(d__1)) * xmaxju[ii - minu];
5865 bid1 = bid0 * xmaxjv[jj - minv] + bid1;
5873 /* ----------------------- Calculate the max error ----------------------*/
5875 bid1 = AdvApp2Var_MathBase::mzsnorm_(ndimen, &vecerr[1]);
5879 *erreur = AdvApp2Var_MathBase::mzsnorm_(&nd, vaux);
5881 /* ------------------------- The end ------------------------------------
5885 AdvApp2Var_SysBase::mgsomsg_("MMA2ER1", 7L);
5890 //=======================================================================
5891 //function : mma2er2_
5893 //=======================================================================
5894 int mma2er2_(integer *ndjacu,
5906 doublereal *epmscut,
5913 /* System generated locals */
5914 integer patjac_dim1, patjac_dim2, patjac_offset, i__1, i__2;
5917 /* Local variables */
5918 static logical ldbg;
5919 static doublereal vaux[2];
5920 static integer i2rdu, i2rdv;
5921 static doublereal errnu, errnv;
5922 static integer ii, nd, jj, nu, nv;
5923 static doublereal bid0, bid1;
5925 /* **********************************************************************
5930 /* Remove coefficients of PATJAC to obtain the minimum degree */
5931 /* by U and V checking the imposed tolerance. */
5935 /* TOUS,AB_SPECIFI:: CARREAU&,CALCUL,&ERREUR */
5937 /* INPUT ARGUMENTS : */
5938 /* ------------------ */
5939 /* NDJACU: Degree by U of table PATJAC. */
5940 /* NDJACV: Degree by V of table PATJAC. */
5941 /* NDIMEN: Dimension of the space. */
5942 /* MINDGU: Limit of index by U of coeff. of PATJAC to be PRESERVED (should be >=0). */
5943 /* MAXDGU: Upper limit of index by U of coeff. of PATJAC to be taken into account. */
5944 /* MINDGV: Limit of index by V of coeff. of PATJAC to be PRESERVED (should be >=0). */
5945 /* MAXDGV: Upper limit of index by V of coeff. of PATJAC to be taken into account. */
5946 /* IORDRU: Order of continuity by U provided by square PATJAC (from -1 to 2) */
5947 /* IORDRV: Order of continuity by U provided by square PATJAC (from -1 to 2) */
5948 /* XMAXJU: Maximum value of Jacobi polynoms of order IORDRU, */
5949 /* from degree 0 to MAXDGU - 2*(IORDU+1) */
5950 /* XMAXJV: Maximum value of Jacobi polynoms of order IORDRV, */
5951 /* from degree 0 to MAXDGV - 2*(IORDV+1) */
5952 /* PATJAC: Table of coeff. of square of approximation with */
5953 /* constraints of order IORDRU by U and IORDRV by V. */
5954 /* EPMSCUT: Tolerance of approximation. */
5955 /* VECERR: Auxiliary vector. */
5956 /* ERREUR: MAX Error commited ALREADY CALCULATED */
5958 /* OUTPUT ARGUMENTS : */
5959 /* ------------------- */
5960 /* ERREUR: MAX Error commited by preserving only coeff of PATJAC */
5961 /* of indices from 0 to NEWDGU by U and from 0 to NEWDGV by V */
5962 /* PLUS the already calculated error. */
5963 /* NEWDGU: Min. Degree by U such as the square of approximation */
5964 /* could check the tolerance. There is always NEWDGU >= MINDGU >= 0. */
5965 /* NEWDGV: Min. Degree by V such as the square of approximation */
5966 /* could check the tolerance. There is always NEWDGV >= MINDGV >= 0. */
5969 /* COMMONS USED : */
5970 /* ---------------- */
5972 /* REFERENCES CALLED : */
5973 /* --------------------- */
5975 /* DESCRIPTION/NOTES/LIMITATIONS : */
5976 /* ----------------------------------- */
5977 /* Table PATJAC is the place of storage of coeff. Cij of the square of */
5978 /* approximation of F(U,V). The indices i and j show the degree */
5979 /* by U and by V of base polynoms. These polynoms have the form: */
5981 /* ((1 - U*U)**(IORDRU+1)).J(i-2*(IORDRU+1)(U), where */
5983 /* polynom J(i-2*(IORDU+1)(U) is the Jacobi polynom of order */
5984 /* IORDRU+1 (the same by V by replacing U u V in the expression above). */
5986 /* The contribution to the error of term Cij when it is */
5987 /* removed from PATJAC is increased by: */
5989 /* DABS(Cij)*XMAXJU(i-2*(IORDRU+1))*XMAXJV(J-2*(IORDRV+1)) where */
5991 /* XMAXJU(i-2*(IORDRU+1) = ((1 - U*U)**(IORDRU+1)).J(i-2*(IORDRU+1)(U),
5993 /* XMAXJV(i-2*(IORDRV+1) = ((1 - V*V)**(IORDRV+1)).J(j-2*(IORDRV+1)(V).
5997 /* **********************************************************************
5999 /* Name of the routine */
6002 /* ----------------------------- Initialisations ------------------------
6005 /* Parameter adjustments */
6007 patjac_dim1 = *ndjacu + 1;
6008 patjac_dim2 = *ndjacv + 1;
6009 patjac_offset = patjac_dim1 * patjac_dim2;
6010 patjac -= patjac_offset;
6013 ldbg = AdvApp2Var_SysBase::mnfndeb_() >= 3;
6015 AdvApp2Var_SysBase::mgenmsg_("MMA2ER2", 7L);
6018 i2rdu = (*iordru + 1) << 1;
6019 i2rdv = (*iordrv + 1) << 1;
6023 /* **********************************************************************
6025 /* -------------------- Cutting of oefficients ------------------------
6027 /* **********************************************************************
6032 /* ------------------- Calculate the increment of max error --------------- */
6033 /* ----- during the removal of coeff. of indices from MINDGU to MAXDGU ------ */
6034 /* ---------------- by U, the degree by V is fixed to NV -----------------
6039 bid0 = xmaxjv[nv - i2rdv];
6041 for (nd = 1; nd <= i__1; ++nd) {
6044 for (ii = i2rdu; ii <= i__2; ++ii) {
6045 bid1 += (d__1 = patjac[ii + (nv + nd * patjac_dim2) *
6046 patjac_dim1], advapp_abs(d__1)) * xmaxju[ii - i2rdu] * bid0;
6053 vecerr[1] = *epmscut * 2;
6055 errnv = AdvApp2Var_MathBase::mzsnorm_(ndimen, &vecerr[1]);
6057 /* ------------------- Calculate the increment of max error --------------- */
6058 /* ----- during the removal of coeff. of indices from MINDGV to MAXDGV ------ */
6059 /* ---------------- by V, the degree by U is fixed to NU -----------------
6064 bid0 = xmaxju[nu - i2rdu];
6066 for (nd = 1; nd <= i__1; ++nd) {
6069 for (jj = i2rdv; jj <= i__2; ++jj) {
6070 bid1 += (d__1 = patjac[nu + (jj + nd * patjac_dim2) *
6071 patjac_dim1], advapp_abs(d__1)) * xmaxjv[jj - i2rdv] * bid0;
6078 vecerr[1] = *epmscut * 2;
6080 errnu = AdvApp2Var_MathBase::mzsnorm_(ndimen, &vecerr[1]);
6082 /* ----------------------- Calculate the max error ----------------------
6088 errnu = AdvApp2Var_MathBase::mzsnorm_(&nd, vaux);
6090 errnv = AdvApp2Var_MathBase::mzsnorm_(&nd, vaux);
6092 if (errnu > errnv) {
6093 if (errnv < *epmscut) {
6100 if (errnu < *epmscut) {
6110 /* -------------------------- Return the degrees -------------------
6114 *newdgu = advapp_max(nu,1);
6115 *newdgv = advapp_max(nv,1);
6117 /* ----------------------------------- The end --------------------------
6121 AdvApp2Var_SysBase::mgsomsg_("MMA2ER2", 7L);
6126 //=======================================================================
6127 //function : mma2fnc_
6129 //=======================================================================
6130 int AdvApp2Var_ApproxF2var::mma2fnc_(integer *ndimen,
6134 const AdvApp2Var_EvaluatorFunc2Var& foncnp,
6158 static integer c__8 = 8;
6160 /* System generated locals */
6161 integer courbe_dim1, courbe_dim2, courbe_offset, somtab_dim1, somtab_dim2,
6162 somtab_offset, diftab_dim1, diftab_dim2, diftab_offset,
6163 contr1_dim1, contr1_dim2, contr1_offset, contr2_dim1, contr2_dim2,
6164 contr2_offset, errmax_dim1, errmax_offset, errmoy_dim1,
6165 errmoy_offset, i__1;
6168 /* Local variables */
6169 static integer ideb;
6170 static doublereal tmil;
6171 static integer ideb1, ibid1, ibid2, ncfja, ndgre, ilong,
6173 static doublereal wrkar[1];
6174 static integer nupil;
6175 static intptr_t iofwr;
6176 static doublereal uvpav[4] /* was [2][2] */;
6177 static integer nd, ii;
6180 static doublereal uv11[4] /* was [2][2] */;
6181 static integer ncb1;
6182 static doublereal eps3;
6183 static integer isz1, isz2, isz3, isz4, isz5;
6184 static intptr_t ipt1, ipt2, ipt3, ipt4, ipt5,iptt, jptt;
6186 /* **********************************************************************
6191 /* Approximation of a limit of non polynomial function F(u,v) */
6192 /* (in the space of dimension NDIMEN) by SEVERAL */
6193 /* polynomial curves, by the method of least squares. The parameter of the function is preserved. */
6197 /* TOUS, AB_SPECIFI :: FONCTION&,EXTREMITE&, APPROXIMATION, &COURBE. */
6199 /* INPUT ARGUMENTS : */
6200 /* ----------------- */
6201 /* NDIMEN: Total Dimension of the space (sum of dimensions */
6202 /* of sub-spaces) */
6203 /* NBSESP: Number of "independent" sub-spaces. */
6204 /* NDIMSE: Table of dimensions of sub-spaces. */
6205 /* UVFONC: Limits of the interval (a,b)x(c,d) of definition of the */
6206 /* function to be approached by U (UVFONC(*,1) contains (a,b)) */
6207 /* and by V (UVFONC(*,2) contains (c,d)). */
6208 /* FONCNP: External function of position on the non polynomial function to be approached. */
6209 /* TCONST: Value of isoparameter of F(u,v) to be discretized. */
6210 /* ISOFAV: Type of chosen iso, = 1, shose that discretization is with u */
6211 /* fixed; = 2, shows that v is fixed. */
6212 /* NBROOT: Nb of points of discretisation of the iso, extremities not included. */
6213 /* ROOTLG: Table of roots of the polynom of Legendre defined on */
6214 /* (-1,1), of degree NBROOT. */
6215 /* IORDRE: Order of constraint at the extremities of the limit */
6216 /* -1 = no constraints, */
6217 /* 0 = constraints of passage to limits (i.e. C0), */
6218 /* 1 = C0 + constraints of 1st derivatives (i.e. C1), */
6219 /* 2 = C1 + constraints of 2nd derivatives (i.e. C2). */
6220 /* IDERIV: Order of derivative of the limit. */
6221 /* NDGJAC: Degree of serial development to be used for calculation in */
6222 /* the Jacobi base. */
6223 /* NBCRMX: Max Nb of curves to be created. */
6224 /* NCFLIM: Max Nb of coeff of the polynomial curve */
6225 /* of approximation (should be above or equal to */
6226 /* 2*IORDRE+2 and below or equal to 50). */
6227 /* EPSAPR: Table of required errors of approximation */
6228 /* sub-space by sub-space. */
6230 /* OUTPUT ARGUMENTS : */
6231 /* ------------------- */
6232 /* NCOEFF: Number of significative coeff of calculated curves. */
6233 /* COURBE: Table of coeff. of calculated polynomial curves. */
6234 /* Should be dimensioned in (NCFLIM,NDIMEN,NBCRMX). */
6235 /* These curves are ALWAYS parametrized in (-1,1). */
6236 /* NBCRBE: Nb of calculated curves. */
6237 /* SOMTAB: For F defined on (-1,1) (otherwise rescale the */
6238 /* parameters), this is the table of sums F(u,vj) + F(u,-vj)
6240 /* if ISOFAV = 1 (and IDERIV=0, otherwise the derivatives */
6241 /* by u of order IDERIV are taken) or sumes F(ui,v) + F(-ui,v) if */
6242 /* ISOFAV = 2 (and IDERIV=0, otherwise the derivatives by */
6243 /* v of order IDERIV are taken). */
6244 /* DIFTAB: For F defined on (-1,1) (otherwise rescale the */
6245 /* parameters), this is the table of sums F(u,vj) - F(u,-vj)
6247 /* if ISOFAV = 1 (and IDERIV=0, otherwise the derivatives */
6248 /* by u of order IDERIV are taken) or sumes F(ui,v) + F(-ui,v) if */
6249 /* ISOFAV = 2 (and IDERIV=0, otherwise the derivatives by */
6250 /* v of order IDERIV are taken). */
6251 /* CONTR1: Contains the coordinates of the left extremity of the iso */
6252 /* and of its derivatives till order IORDRE */
6253 /* CONTR2: Contains the coordinates of the right extremity of the iso */
6254 /* and of its derivatives till order IORDRE */
6255 /* TABDEC: Table of NBCRBE+1 parameters of cut of UVFONC(1:2,1)
6257 /* if ISOFAV=2, or of UVFONC(1:2,2) if ISOFAV=1. */
6258 /* ERRMAX: Table of MAX errors (sub-space by sub-space) */
6259 /* committed in the approximation of FONCNP by NBCRBE curves. */
6260 /* ERRMOY: Table of AVERAGE errors (sub-space by sub-space) */
6261 /* committed in the approximation of FONCNP by NBCRBE curves.
6262 /* IERCOD: Error code: */
6263 /* -1 = ERRMAX > EPSAPR for at least one sub-space. */
6264 /* (the resulting curves of at least mathematic degree NCFLIM-1 */
6265 /* are calculated). */
6266 /* 0 = Everything is ok. */
6267 /* 1 = Pb of incoherence of inputs. */
6268 /* 10 = Pb of calculation of the interpolation of constraints. */
6269 /* 13 = Pb in the dynamic allocation. */
6270 /* 33 = Pb in the data recuperation from block data */
6271 /* of coeff. of integration by GAUSS method. */
6272 /* >100 Pb in the evaluation of FONCNP, the returned error code */
6273 /* is equal to the error code of FONCNP + 100. */
6275 /* COMMONS USED : */
6276 /* ---------------- */
6278 /* REFERENCES CALLED : */
6279 /* ----------------------- */
6281 /* DESCRIPTION/NOTES/LIMITATIONS : */
6282 /* ----------------------------------- */
6283 /* --> The approximation part is done in the space of dimension */
6284 /* NDIMEN (the sum of dimensions of sub-spaces). For example : */
6285 /* If NBSESP=2 and NDIMSE(1)=3, NDIMSE(2)=2, there is smoothing with */
6286 /* NDIMEN=5. The result (in COURBE(NDIMEN,NCOEFF,i) ), will be */
6287 /* composed of the result of smoothing of 3D function in */
6288 /* COURBE(1:3,1:NCOEFF,i) and of smoothing of 2D function in */
6289 /* COURBE(4:5,1:NCOEFF,i). */
6291 /* --> Routine FONCNP should be declared EXTERNAL in the program */
6292 /* calling MMA2FNC. */
6294 /* --> Function FONCNP, declared externally, should be declared */
6295 /* IMPERATIVELY in form : */
6296 /* SUBROUTINE FONCNP(NDIMEN,UINTFN,VINTFN,IIUOUV,TCONST,NBPTAB */
6297 /* ,TTABLE,IDERIU,IDERIV,IERCOD) */
6298 /* where the input arguments are : */
6299 /* - NDIMEN is integer defined as the sum of dimensions of */
6300 /* sub-spaces (i.e. total dimension of the problem). */
6301 /* - UINTFN(2) is a table of 2 reals containing the interval */
6302 /* by u where the function to be approximated is defined */
6303 /* (so it is equal to UIFONC). */
6304 /* - VINTFN(2) is a table of 2 reals containing the interval */
6305 /* by v where the function to be approximated is defined */
6306 /* (so it is equal to VIFONC). */
6307 /* - IIUOUV, shows that the points to be calculated have a constant U */
6308 /* (IIUOUV=1) or a constant V (IIUOUV=2). */
6309 /* - TCONST, real, value of the fixed discretisation parameter. Takes values */
6310 /* in (UINTFN(1),UINTFN(2)) if IIUOUV=1, */
6311 /* or in (VINTFN(1),VINTFN(2)) if IIUOUV=2. */
6312 /* - NBPTAB, the nb of point of discretisation following the free variable */
6313 /* : V if IIUOUV=1 or U if IIUOUV = 2. */
6314 /* - TTABLE, Table of NBPTAB parametres of discretisation. . */
6315 /* - IDERIU, integer, takes values between 0 (position) */
6316 /* and IORDREU (partial derivative of the function by u */
6317 /* of order IORDREU if IORDREU > 0). */
6318 /* - IDERIV, integer, takes values between 0 (position) */
6319 /* and IORDREV (partial derivative of the function by v */
6320 /* of order IORDREV if IORDREV > 0). */
6321 /* and the output arguments are : */
6322 /* - FPNTAB(NDIMEN,NBPTAB) contains, at output, the table of */
6323 /* NBPTAB points calculated in FONCNP. */
6324 /* - IERCOD is, at output the error code of FONCNP. This code */
6325 /* (integer) should be strictly positive if there is a problem. */
6327 /* The input arguments SHOULD NOT BE modified under FONCNP.
6330 /* --> If IERCOD=-1, the required precision can't be reached (ERRMAX */
6331 /* is above EPSAPR on at least one sub-space), but
6333 /* one gives the best possible result for NCFLIM and EPSAPR */
6334 /* chosen by the user. In this case (and for IERCOD=0), there is a solution. */
6337 /* **********************************************************************
6339 /* Name of the routine */
6341 /* Parameter adjustments */
6346 errmoy_dim1 = *nbsesp;
6347 errmoy_offset = errmoy_dim1 + 1;
6348 errmoy -= errmoy_offset;
6349 errmax_dim1 = *nbsesp;
6350 errmax_offset = errmax_dim1 + 1;
6351 errmax -= errmax_offset;
6352 contr2_dim1 = *ndimen;
6353 contr2_dim2 = *iordre + 2;
6354 contr2_offset = contr2_dim1 * (contr2_dim2 + 1) + 1;
6355 contr2 -= contr2_offset;
6356 contr1_dim1 = *ndimen;
6357 contr1_dim2 = *iordre + 2;
6358 contr1_offset = contr1_dim1 * (contr1_dim2 + 1) + 1;
6359 contr1 -= contr1_offset;
6360 diftab_dim1 = *nbroot / 2 + 1;
6361 diftab_dim2 = *ndimen;
6362 diftab_offset = diftab_dim1 * (diftab_dim2 + 1);
6363 diftab -= diftab_offset;
6364 somtab_dim1 = *nbroot / 2 + 1;
6365 somtab_dim2 = *ndimen;
6366 somtab_offset = somtab_dim1 * (somtab_dim2 + 1);
6367 somtab -= somtab_offset;
6369 courbe_dim1 = *ncflim;
6370 courbe_dim2 = *ndimen;
6371 courbe_offset = courbe_dim1 * (courbe_dim2 + 1) + 1;
6372 courbe -= courbe_offset;
6375 ibb = AdvApp2Var_SysBase::mnfndeb_();
6377 AdvApp2Var_SysBase::mgenmsg_("MMA2FNC", 7L);
6382 /* ---------------- Set to zero the coefficients of CURVE --------------
6385 ilong = *ndimen * *ncflim * *nbcrmx;
6386 AdvApp2Var_SysBase::mvriraz_(&ilong, &courbe[courbe_offset]);
6388 /* **********************************************************************
6390 /* -------------------------- Checking of entries ------------------
6392 /* **********************************************************************
6395 AdvApp2Var_MathBase::mmveps3_(&eps3);
6396 if ((d__1 = uvfonc[4] - uvfonc[3], advapp_abs(d__1)) < eps3) {
6399 if ((d__1 = uvfonc[6] - uvfonc[5], advapp_abs(d__1)) < eps3) {
6408 /* ********************************************************************** */
6409 /* ------------- Preparation of parameters of discretisation ----------- */
6410 /* **********************************************************************
6413 /* -- Allocation of a table of parameters and points of discretisation -- */
6414 /* --> For the parameters of discretisation. */
6416 /* --> For the points of discretisation in MMA1FDI and MMA1CDI and
6418 /* the auxiliary curve for MMAPCMP */
6419 ibid1 = *ndimen * (*nbroot + 2);
6420 ibid2 = ((*iordre + 1) << 1) * *nbroot;
6421 isz2 = advapp_max(ibid1,ibid2);
6422 ibid1 = (((*ncflim - 1) / 2 + 1) << 1) * *ndimen;
6423 isz2 = advapp_max(ibid1,isz2);
6424 /* --> To return the polynoms of hermit. */
6425 isz3 = ((*iordre + 1) << 2) * (*iordre + 1);
6426 /* --> For the Gauss coeff. of integration. */
6427 isz4 = (*nbroot / 2 + 1) * (*ndgjac + 1 - ((*iordre + 1) << 1));
6428 /* --> For the coeff of the curve in the base of Jacobi */
6429 isz5 = (*ndgjac + 1) * *ndimen;
6431 ndwrk = isz1 + isz2 + isz3 + isz4 + isz5;
6432 AdvApp2Var_SysBase::mcrrqst_(&c__8, &ndwrk, wrkar, &iofwr, &ier);
6435 /* --> For the parameters of discretisation (NBROOT+2 extremities). */
6437 /* --> For the points of discretisation FPNTAB(NDIMEN,NBROOT+2), */
6438 /* FPNTAB(NBROOT,2*(IORDRE+1)) and for WRKAR of MMAPCMP. */
6440 /* --> For the polynoms of Hermit */
6442 /* --> For the Gauss coeff of integration. */
6444 /* --> For the curve in Jacobi. */
6447 /* ------------------ Initialisation of management of cuts ---------
6451 uvpav[0] = uvfonc[3];
6452 uvpav[1] = uvfonc[4];
6453 tabdec[0] = uvfonc[5];
6454 tabdec[1] = uvfonc[6];
6455 } else if (*isofav == 2) {
6456 tabdec[0] = uvfonc[3];
6457 tabdec[1] = uvfonc[4];
6458 uvpav[2] = uvfonc[5];
6459 uvpav[3] = uvfonc[6];
6467 /* **********************************************************************
6469 /* APPROXIMATION WITH CUTS */
6470 /* **********************************************************************
6474 /* --> When the top is reached, this is the end ! */
6475 if (nupil - *nbcrbe == 0) {
6480 uvpav[2] = tabdec[*nbcrbe];
6481 uvpav[3] = tabdec[*nbcrbe + 1];
6482 } else if (*isofav == 2) {
6483 uvpav[0] = tabdec[*nbcrbe];
6484 uvpav[1] = tabdec[*nbcrbe + 1];
6489 /* -------------------- Normalization of parameters -------------------- */
6491 mma1nop_(nbroot, &rootlg[1], uvpav, isofav, &wrkar[ipt1], &ier);
6496 /* -------------------- Discretisation of FONCNP ------------------------ */
6498 mma1fdi_(ndimen, uvpav, foncnp, isofav, tconst, nbroot, &wrkar[ipt1],
6499 iordre, ideriv, &wrkar[ipt2], &somtab[(ncb1 * somtab_dim2 + 1) *
6500 somtab_dim1], &diftab[(ncb1 * diftab_dim2 + 1) * diftab_dim1], &
6501 contr1[(ncb1 * contr1_dim2 + 1) * contr1_dim1 + 1], &contr2[(ncb1
6502 * contr2_dim2 + 1) * contr2_dim1 + 1], iercod);
6507 /* -----------Cut the discretisation of constraints ------------*/
6510 mma1cdi_(ndimen, nbroot, &rootlg[1], iordre, &contr1[(ncb1 *
6511 contr1_dim2 + 1) * contr1_dim1 + 1], &contr2[(ncb1 *
6512 contr2_dim2 + 1) * contr2_dim1 + 1], &somtab[(ncb1 *
6513 somtab_dim2 + 1) * somtab_dim1], &diftab[(ncb1 * diftab_dim2
6514 + 1) * diftab_dim1], &wrkar[ipt2], &wrkar[ipt3], &ier);
6520 /* **********************************************************************
6522 /* -------------------- Calculate the curve of approximation -------------
6524 /* **********************************************************************
6527 mma1jak_(ndimen, nbroot, iordre, ndgjac, &somtab[(ncb1 * somtab_dim2 + 1)
6528 * somtab_dim1], &diftab[(ncb1 * diftab_dim2 + 1) * diftab_dim1], &
6529 wrkar[ipt4], &wrkar[ipt5], &ier);
6534 /* **********************************************************************
6536 /* ---------------- Add polynom of interpolation -------------------
6538 /* **********************************************************************
6542 mma1cnt_(ndimen, iordre, &contr1[(ncb1 * contr1_dim2 + 1) *
6543 contr1_dim1 + 1], &contr2[(ncb1 * contr2_dim2 + 1) *
6544 contr2_dim1 + 1], &wrkar[ipt3], ndgjac, &wrkar[ipt5]);
6547 /* **********************************************************************
6549 /* --------------- Calculate Max and Average error ----------------------
6551 /* **********************************************************************
6554 mma1fer_(ndimen, nbsesp, &ndimse[1], iordre, ndgjac, &wrkar[ipt5], ncflim,
6555 &epsapr[1], &wrkar[ipt2], &errmax[ncb1 * errmax_dim1 + 1], &
6556 errmoy[ncb1 * errmoy_dim1 + 1], &ncoeff[ncb1], &ier);
6561 if (ier == 0 || (ier == -1 && nupil == *nbcrmx)) {
6563 /* ******************************************************************
6565 /* ----------------------- Compression du resultat ------------------
6567 /* ******************************************************************
6573 ncfja = *ndgjac + 1;
6574 /* -> Compression of result in WRKAR(IPT2) */
6577 AdvApp2Var_MathBase::mmapcmp_(ndimen,
6578 &ncfja, &ncoeff[ncb1], &wrkar[ipt5], &wrkar[ipt2]);
6580 AdvApp2Var_MathBase::mmapcmp_((integer*)ndimen,
6586 ilong = *ndimen * *ncflim;
6587 AdvApp2Var_SysBase::mvriraz_(&ilong, &wrkar[ipt5]);
6588 /* -> Passage to canonic base (-1,1) (result in WRKAR(IPT5)).
6590 ndgre = ncoeff[ncb1] - 1;
6592 for (nd = 1; nd <= i__1; ++nd) {
6593 iptt = ipt2 + ((nd - 1) << 1) * (ndgre / 2 + 1);
6594 jptt = ipt5 + (nd - 1) * ncoeff[ncb1];
6595 AdvApp2Var_MathBase::mmjacan_(iordre, &ndgre, &wrkar[iptt], &wrkar[jptt]);
6599 /* -> Store the calculated curve */
6601 AdvApp2Var_MathBase::mmfmca8_(&ncoeff[ncb1], ndimen, &ibid1, ncflim, ndimen, &ibid1, &
6602 wrkar[ipt5], &courbe[(ncb1 * courbe_dim2 + 1) * courbe_dim1 +
6605 /* -> Before normalization of constraints on (-1,1), recalculate */
6606 /* the true constraints. */
6608 for (ii = 0; ii <= i__1; ++ii) {
6609 mma1noc_(uv11, ndimen, &ii, &contr1[(ii + 1 + ncb1 * contr1_dim2)
6610 * contr1_dim1 + 1], uvpav, isofav, ideriv, &contr1[(ii +
6611 1 + ncb1 * contr1_dim2) * contr1_dim1 + 1]);
6612 mma1noc_(uv11, ndimen, &ii, &contr2[(ii + 1 + ncb1 * contr2_dim2)
6613 * contr2_dim1 + 1], uvpav, isofav, ideriv, &contr2[(ii +
6614 1 + ncb1 * contr2_dim2) * contr2_dim1 + 1]);
6618 ibid1 = (*nbroot / 2 + 1) * *ndimen;
6619 mma1noc_(uv11, &ibid1, &ii, &somtab[(ncb1 * somtab_dim2 + 1) *
6620 somtab_dim1], uvpav, isofav, ideriv, &somtab[(ncb1 *
6621 somtab_dim2 + 1) * somtab_dim1]);
6622 mma1noc_(uv11, &ibid1, &ii, &diftab[(ncb1 * diftab_dim2 + 1) *
6623 diftab_dim1], uvpav, isofav, ideriv, &diftab[(ncb1 *
6624 diftab_dim2 + 1) * diftab_dim1]);
6627 for (nd = 1; nd <= i__1; ++nd) {
6628 mma1noc_(uv11, &ncoeff[ncb1], &ii, &courbe[(nd + ncb1 *
6629 courbe_dim2) * courbe_dim1 + 1], uvpav, isofav, ideriv, &
6630 courbe[(nd + ncb1 * courbe_dim2) * courbe_dim1 + 1]);
6634 /* -> Update the nb of already created curves */
6637 /* -> ...otherwise try to cut the current interval in 2... */
6639 tmil = (tabdec[*nbcrbe + 1] + tabdec[*nbcrbe]) / 2.;
6642 ilong = (nupil - *nbcrbe) << 3;
6643 AdvApp2Var_SysBase::mcrfill_(&ilong, &tabdec[ideb],&tabdec[ideb1]);
6644 tabdec[ideb] = tmil;
6648 /* ---------- Make approximation of the rest -----------
6653 /* --------------------- Return code of error -----------------
6655 /* --> Pb with dynamic allocation */
6659 /* --> Inputs incoherent. */
6664 /* -------------------------- Dynamic desallocation -------------------
6669 AdvApp2Var_SysBase::mcrdelt_(&c__8, &ndwrk, wrkar, &iofwr, &ier);
6676 /* ------------------------------ The end -------------------------------
6681 AdvApp2Var_SysBase::maermsg_("MMA2FNC", iercod, 7L);
6684 AdvApp2Var_SysBase::mgsomsg_("MMA2FNC", 7L);
6689 //=======================================================================
6690 //function : mma2fx6_
6692 //=======================================================================
6693 int AdvApp2Var_ApproxF2var::mma2fx6_(integer *ncfmxu,
6710 /* System generated locals */
6711 integer epsfro_dim1, epsfro_offset, patcan_dim1, patcan_dim2, patcan_dim3,
6712 patcan_dim4, patcan_offset, errmax_dim1, errmax_dim2,
6713 errmax_offset, ncoefu_dim1, ncoefu_offset, ncoefv_dim1,
6714 ncoefv_offset, i__1, i__2, i__3, i__4, i__5;
6715 doublereal d__1, d__2;
6717 /* Local variables */
6718 static integer idim, ncfu, ncfv, id, ii, nd, jj, ku, kv, ns, ibb;
6719 static doublereal bid;
6720 static doublereal tol;
6722 /* **********************************************************************
6727 /* Reduction of degree when the squares are the squares of constraints. */
6731 /* TOUS,AB_SPECIFI::CARREAU&,REDUCTION,&CARREAU */
6733 /* INPUT ARGUMENTS : */
6734 /* ------------------ */
6735 /* NCFMXU: Max Nb of coeff by u of solution P(u,v) (table */
6736 /* PATCAN). This argument serves only to declare the size of this table. */
6737 /* NCFMXV: Max Nb of coeff by v of solution P(u,v) (table */
6738 /* PATCAN). This argument serves only to declare the size of this table. */
6739 /* NDIMEN: Total dimension of the space where the processed function */
6740 /* takes its values.(sum of dimensions of sub-spaces) */
6741 /* NBSESP: Nb of independent sub-spaces where the errors are measured. */
6742 /* NDIMSE: Table of dimensions of NBSESP sub-spaces. */
6743 /* NBUPAT: Nb of square solution by u. */
6744 /* NBVPAT: Nb of square solution by v. */
6745 /* IORDRU: Order of constraint imposed at the extremities of iso-V */
6746 /* = 0, the extremities of iso-V are calculated */
6747 /* = 1, additionally the 1st derivative in the direction of iso-V is calculated */
6748 /* = 2, additionally the 2nd derivative in the direction of iso-V is calculated */
6749 /* IORDRV: Ordre de contrainte impose aux extremites de l'iso-U */
6750 /* = 0, on calcule les extremites de l'iso-U. */
6751 /* = 1, additionally the 1st derivative in the direction of iso-U is calculated */
6752 /* = 2, additionally the 2nd derivative in the direction of iso-U is calculated */
6753 /* EPSAPR: Table of imposed precisions, sub-space by sub-space. */
6754 /* EPSFRO: Table of imposed precisions, sub-space by sub-space on the limits of squares. */
6755 /* PATCAN: Table of coeff. in the canonic base of squares P(u,v) calculated for (u,v) in (-1,1). */
6756 /* ERRMAX: Table of MAX errors (sub-space by sub-space) */
6757 /* committed in the approximation of F(u,v) by P(u,v). */
6758 /* NCOEFU: Table of Nb of significative coeffs. by u of calculated squares. */
6759 /* NCOEFV: Table of Nb of significative coeffs. by v of calculated squares. */
6761 /* OUTPUT ARGUMENTS : */
6762 /* ------------------- */
6763 /* NCOEFU: Table of Nb of significative coeffs. by u of calculated squares. */
6764 /* NCOEFV: Table of Nb of significative coeffs. by v of calculated squares. */
6766 /* COMMONS USED : */
6767 /* ---------------- */
6769 /* REFERENCES CALLED : */
6770 /* --------------------- */
6772 /* DESCRIPTION/NOTES/LIMITATIONS : */
6773 /* ------------------------------- */
6775 /* **********************************************************************
6778 /* Name of the routine */
6781 /* Parameter adjustments */
6782 epsfro_dim1 = *nbsesp;
6783 epsfro_offset = epsfro_dim1 * 5 + 1;
6784 epsfro -= epsfro_offset;
6787 ncoefv_dim1 = *nbupat;
6788 ncoefv_offset = ncoefv_dim1 + 1;
6789 ncoefv -= ncoefv_offset;
6790 ncoefu_dim1 = *nbupat;
6791 ncoefu_offset = ncoefu_dim1 + 1;
6792 ncoefu -= ncoefu_offset;
6793 errmax_dim1 = *nbsesp;
6794 errmax_dim2 = *nbupat;
6795 errmax_offset = errmax_dim1 * (errmax_dim2 + 1) + 1;
6796 errmax -= errmax_offset;
6797 patcan_dim1 = *ncfmxu;
6798 patcan_dim2 = *ncfmxv;
6799 patcan_dim3 = *ndimen;
6800 patcan_dim4 = *nbupat;
6801 patcan_offset = patcan_dim1 * (patcan_dim2 * (patcan_dim3 * (patcan_dim4
6803 patcan -= patcan_offset;
6806 ibb = AdvApp2Var_SysBase::mnfndeb_();
6808 AdvApp2Var_SysBase::mgenmsg_("MMA2FX6", 7L);
6813 for (jj = 1; jj <= i__1; ++jj) {
6815 for (ii = 1; ii <= i__2; ++ii) {
6816 ncfu = ncoefu[ii + jj * ncoefu_dim1];
6817 ncfv = ncoefv[ii + jj * ncoefv_dim1];
6819 /* ********************************************************************** */
6820 /* -------------------- Reduction of degree by U ------------------------- */
6821 /* ********************************************************************** */
6824 if (ncfu <= (*iordru + 1) << 1 && ncfu > 2) {
6828 for (ns = 1; ns <= i__3; ++ns) {
6831 d__1 = tol, d__2 = epsfro[ns + epsfro_dim1 * 9];
6832 tol = advapp_min(d__1,d__2);
6834 d__1 = tol, d__2 = epsfro[ns + epsfro_dim1 * 10];
6835 tol = advapp_min(d__1,d__2);
6837 d__1 = tol, d__2 = epsfro[ns + epsfro_dim1 * 11];
6838 tol = advapp_min(d__1,d__2);
6840 d__1 = tol, d__2 = epsfro[ns + epsfro_dim1 * 12];
6841 tol = advapp_min(d__1,d__2);
6842 if (ii == 1 || ii == *nbupat || jj == 1 || jj == *nbvpat)
6845 d__1 = tol, d__2 = epsfro[ns + epsfro_dim1 * 5];
6846 tol = advapp_min(d__1,d__2);
6848 d__1 = tol, d__2 = epsfro[ns + epsfro_dim1 * 6];
6849 tol = advapp_min(d__1,d__2);
6851 d__1 = tol, d__2 = epsfro[ns + epsfro_dim1 * 7];
6852 tol = advapp_min(d__1,d__2);
6854 d__1 = tol, d__2 = epsfro[ns + (epsfro_dim1 << 3)];
6855 tol = advapp_min(d__1,d__2);
6860 for (nd = 1; nd <= i__4; ++nd) {
6863 for (kv = 1; kv <= i__5; ++kv) {
6864 bid += (d__1 = patcan[ncfu + (kv + (id + (ii + jj
6865 * patcan_dim4) * patcan_dim3) *
6866 patcan_dim2) * patcan_dim1], advapp_abs(d__1));
6872 if (bid > tol * 1e-6 || bid > errmax[ns + (ii + jj *
6873 errmax_dim2) * errmax_dim1]) {
6884 /* ********************************************************************** */
6885 /* -------------------- Reduction of degree by V ------------------------- */
6886 /* ********************************************************************** */
6889 if (ncfv <= (*iordrv + 1) << 1 && ncfv > 2) {
6893 for (ns = 1; ns <= i__3; ++ns) {
6896 d__1 = tol, d__2 = epsfro[ns + epsfro_dim1 * 9];
6897 tol = advapp_min(d__1,d__2);
6899 d__1 = tol, d__2 = epsfro[ns + epsfro_dim1 * 10];
6900 tol = advapp_min(d__1,d__2);
6902 d__1 = tol, d__2 = epsfro[ns + epsfro_dim1 * 11];
6903 tol = advapp_min(d__1,d__2);
6905 d__1 = tol, d__2 = epsfro[ns + epsfro_dim1 * 12];
6906 tol = advapp_min(d__1,d__2);
6907 if (ii == 1 || ii == *nbupat || jj == 1 || jj == *nbvpat)
6910 d__1 = tol, d__2 = epsfro[ns + epsfro_dim1 * 5];
6911 tol = advapp_min(d__1,d__2);
6913 d__1 = tol, d__2 = epsfro[ns + epsfro_dim1 * 6];
6914 tol = advapp_min(d__1,d__2);
6916 d__1 = tol, d__2 = epsfro[ns + epsfro_dim1 * 7];
6917 tol = advapp_min(d__1,d__2);
6919 d__1 = tol, d__2 = epsfro[ns + (epsfro_dim1 << 3)];
6920 tol = advapp_min(d__1,d__2);
6925 for (nd = 1; nd <= i__4; ++nd) {
6928 for (ku = 1; ku <= i__5; ++ku) {
6929 bid += (d__1 = patcan[ku + (ncfv + (id + (ii + jj
6930 * patcan_dim4) * patcan_dim3) *
6931 patcan_dim2) * patcan_dim1], advapp_abs(d__1));
6937 if (bid > tol * 1e-6 || bid > errmax[ns + (ii + jj *
6938 errmax_dim2) * errmax_dim1]) {
6949 /* --- Return the nbs of coeff. and pass to the next square --- */
6952 ncoefu[ii + jj * ncoefu_dim1] = advapp_max(ncfu,2);
6953 ncoefv[ii + jj * ncoefv_dim1] = advapp_max(ncfv,2);
6959 /* ------------------------------ The End -------------------------------
6963 AdvApp2Var_SysBase::mgsomsg_("MMA2FX6", 7L);
6969 //=======================================================================
6970 //function : mma2jmx_
6972 //=======================================================================
6973 int AdvApp2Var_ApproxF2var::mma2jmx_(integer *ndgjac,
6977 /* Initialized data */
6979 static doublereal xmax2[57] = { .9682458365518542212948163499456,
6980 .986013297183269340427888048593603,
6981 1.07810420343739860362585159028115,
6982 1.17325804490920057010925920756025,
6983 1.26476561266905634732910520370741,
6984 1.35169950227289626684434056681946,
6985 1.43424378958284137759129885012494,
6986 1.51281316274895465689402798226634,
6987 1.5878364329591908800533936587012,
6988 1.65970112228228167018443636171226,
6989 1.72874345388622461848433443013543,
6990 1.7952515611463877544077632304216,
6991 1.85947199025328260370244491818047,
6992 1.92161634324190018916351663207101,
6993 1.98186713586472025397859895825157,
6994 2.04038269834980146276967984252188,
6995 2.09730119173852573441223706382076,
6996 2.15274387655763462685970799663412,
6997 2.20681777186342079455059961912859,
6998 2.25961782459354604684402726624239,
6999 2.31122868752403808176824020121524,
7000 2.36172618435386566570998793688131,
7001 2.41117852396114589446497298177554,
7002 2.45964731268663657873849811095449,
7003 2.50718840313973523778244737914028,
7004 2.55385260994795361951813645784034,
7005 2.59968631659221867834697883938297,
7006 2.64473199258285846332860663371298,
7007 2.68902863641518586789566216064557,
7008 2.73261215675199397407027673053895,
7009 2.77551570192374483822124304745691,
7010 2.8177699459714315371037628127545,
7011 2.85940333797200948896046563785957,
7012 2.90044232019793636101516293333324,
7013 2.94091151970640874812265419871976,
7014 2.98083391718088702956696303389061,
7015 3.02023099621926980436221568258656,
7016 3.05912287574998661724731962377847,
7017 3.09752842783622025614245706196447,
7018 3.13546538278134559341444834866301,
7019 3.17295042316122606504398054547289,
7020 3.2099992681699613513775259670214,
7021 3.24662674946606137764916854570219,
7022 3.28284687953866689817670991319787,
7023 3.31867291347259485044591136879087,
7024 3.35411740487202127264475726990106,
7025 3.38919225660177218727305224515862,
7026 3.42390876691942143189170489271753,
7027 3.45827767149820230182596660024454,
7028 3.49230918177808483937957161007792,
7029 3.5260130200285724149540352829756,
7030 3.55939845146044235497103883695448,
7031 3.59247431368364585025958062194665,
7032 3.62524904377393592090180712976368,
7033 3.65773070318071087226169680450936,
7034 3.68992700068237648299565823810245,
7035 3.72184531357268220291630708234186 };
7036 static doublereal xmax4[55] = { 1.1092649593311780079813740546678,
7037 1.05299572648705464724876659688996,
7038 1.0949715351434178709281698645813,
7039 1.15078388379719068145021100764647,
7040 1.2094863084718701596278219811869,
7041 1.26806623151369531323304177532868,
7042 1.32549784426476978866302826176202,
7043 1.38142537365039019558329304432581,
7044 1.43575531950773585146867625840552,
7045 1.48850442653629641402403231015299,
7046 1.53973611681876234549146350844736,
7047 1.58953193485272191557448229046492,
7048 1.63797820416306624705258190017418,
7049 1.68515974143594899185621942934906,
7050 1.73115699602477936547107755854868,
7051 1.77604489805513552087086912113251,
7052 1.81989256661534438347398400420601,
7053 1.86276344480103110090865609776681,
7054 1.90471563564740808542244678597105,
7055 1.94580231994751044968731427898046,
7056 1.98607219357764450634552790950067,
7057 2.02556989246317857340333585562678,
7058 2.06433638992049685189059517340452,
7059 2.10240936014742726236706004607473,
7060 2.13982350649113222745523925190532,
7061 2.17661085564771614285379929798896,
7062 2.21280102016879766322589373557048,
7063 2.2484214321456956597803794333791,
7064 2.28349755104077956674135810027654,
7065 2.31805304852593774867640120860446,
7066 2.35210997297725685169643559615022,
7067 2.38568889602346315560143377261814,
7068 2.41880904328694215730192284109322,
7069 2.45148841120796359750021227795539,
7070 2.48374387161372199992570528025315,
7071 2.5155912654873773953959098501893,
7072 2.54704548720896557684101746505398,
7073 2.57812056037881628390134077704127,
7074 2.60882970619319538196517982945269,
7075 2.63918540521920497868347679257107,
7076 2.66919945330942891495458446613851,
7077 2.69888301230439621709803756505788,
7078 2.72824665609081486737132853370048,
7079 2.75730041251405791603760003778285,
7080 2.78605380158311346185098508516203,
7081 2.81451587035387403267676338931454,
7082 2.84269522483114290814009184272637,
7083 2.87060005919012917988363332454033,
7084 2.89823818258367657739520912946934,
7085 2.92561704377132528239806135133273,
7086 2.95274375377994262301217318010209,
7087 2.97962510678256471794289060402033,
7088 3.00626759936182712291041810228171,
7089 3.03267744830655121818899164295959,
7090 3.05886060707437081434964933864149 };
7091 static doublereal xmax6[53] = { 1.21091229812484768570102219548814,
7092 1.11626917091567929907256116528817,
7093 1.1327140810290884106278510474203,
7094 1.1679452722668028753522098022171,
7095 1.20910611986279066645602153641334,
7096 1.25228283758701572089625983127043,
7097 1.29591971597287895911380446311508,
7098 1.3393138157481884258308028584917,
7099 1.3821288728999671920677617491385,
7100 1.42420414683357356104823573391816,
7101 1.46546895108549501306970087318319,
7102 1.50590085198398789708599726315869,
7103 1.54550385142820987194251585145013,
7104 1.58429644271680300005206185490937,
7105 1.62230484071440103826322971668038,
7106 1.65955905239130512405565733793667,
7107 1.69609056468292429853775667485212,
7108 1.73193098017228915881592458573809,
7109 1.7671112206990325429863426635397,
7110 1.80166107681586964987277458875667,
7111 1.83560897003644959204940535551721,
7112 1.86898184653271388435058371983316,
7113 1.90180515174518670797686768515502,
7114 1.93410285411785808749237200054739,
7115 1.96589749778987993293150856865539,
7116 1.99721027139062501070081653790635,
7117 2.02806108474738744005306947877164,
7118 2.05846864831762572089033752595401,
7119 2.08845055210580131460156962214748,
7120 2.11802334209486194329576724042253,
7121 2.14720259305166593214642386780469,
7122 2.17600297710595096918495785742803,
7123 2.20443832785205516555772788192013,
7124 2.2325216999457379530416998244706,
7125 2.2602654243075083168599953074345,
7126 2.28768115912702794202525264301585,
7127 2.3147799369092684021274946755348,
7128 2.34157220782483457076721300512406,
7129 2.36806787963276257263034969490066,
7130 2.39427635443992520016789041085844,
7131 2.42020656255081863955040620243062,
7132 2.44586699364757383088888037359254,
7133 2.47126572552427660024678584642791,
7134 2.49641045058324178349347438430311,
7135 2.52130850028451113942299097584818,
7136 2.54596686772399937214920135190177,
7137 2.5703922285006754089328998222275,
7138 2.59459096001908861492582631591134,
7139 2.61856915936049852435394597597773,
7140 2.64233265984385295286445444361827,
7141 2.66588704638685848486056711408168,
7142 2.68923766976735295746679957665724,
7143 2.71238965987606292679677228666411 };
7145 /* System generated locals */
7148 /* Local variables */
7149 static logical ldbg;
7150 static integer numax, ii;
7151 static doublereal bid;
7154 /* **********************************************************************
7159 /* Calculate the max of Jacobo polynoms multiplied by the weight on */
7160 /* (-1,1) for order 0,4,6 or Legendre. */
7164 /* LEGENDRE,APPROXIMATION,ERREUR. */
7166 /* INPUT ARGUMENTS : */
7167 /* ------------------ */
7168 /* NDGJAC: Nb of Jacobi coeff. of approximation. */
7169 /* IORDRE: Order of continuity (from -1 to 2) */
7171 /* OUTPUT ARGUMENTS : */
7172 /* ------------------- */
7173 /* XJACMX: Table of maximums of Jacobi polynoms. */
7175 /* COMMONS USED : */
7176 /* ---------------- */
7178 /* REFERENCES CALLED : */
7179 /* --------------------- */
7181 /* DESCRIPTION/NOTES/LIMITATIONS : */
7182 /* ----------------------------------- */
7185 /* ***********************************************************************
7187 /* Name of the routine */
7188 /* ----------------------------- Initialisations ------------------------
7191 ldbg = AdvApp2Var_SysBase::mnfndeb_() >= 3;
7193 AdvApp2Var_SysBase::mgenmsg_("MMA2JMX", 7L);
7196 numax = *ndgjac - ((*iordre + 1) << 1);
7197 if (*iordre == -1) {
7199 for (ii = 0; ii <= i__1; ++ii) {
7200 bid = (ii * 2. + 1.) / 2.;
7201 xjacmx[ii] = sqrt(bid);
7204 } else if (*iordre == 0) {
7206 for (ii = 0; ii <= i__1; ++ii) {
7207 xjacmx[ii] = xmax2[ii];
7210 } else if (*iordre == 1) {
7212 for (ii = 0; ii <= i__1; ++ii) {
7213 xjacmx[ii] = xmax4[ii];
7216 } else if (*iordre == 2) {
7218 for (ii = 0; ii <= i__1; ++ii) {
7219 xjacmx[ii] = xmax6[ii];
7224 /* ------------------------- The end ------------------------------------
7228 AdvApp2Var_SysBase::mgsomsg_("MMA2JMX", 7L);
7233 //=======================================================================
7234 //function : mma2moy_
7236 //=======================================================================
7237 int mma2moy_(integer *ndgumx,
7249 /* System generated locals */
7250 integer patjac_dim1, patjac_dim2, patjac_offset, i__1, i__2, i__3;
7252 /* Local variables */
7253 static logical ldbg;
7254 static integer minu, minv, idebu, idebv, ii, nd, jj;
7255 static doublereal bid0, bid1;
7258 /* **********************************************************************
7263 /* Calculate the average approximation error made when only */
7264 /* the coefficients of PATJAC of degree between */
7265 /* 2*(IORDRU+1) and MINDGU by U and 2*(IORDRV+1) and MINDGV by V are preserved. */
7269 /* LEGENDRE,APPROXIMATION, AVERAGE ERROR */
7271 /* INPUT ARGUMENTS : */
7272 /* ------------------ */
7273 /* NDGUMX: Dimension by U of table PATJAC. */
7274 /* NDGVMX: Dimension by V of table PATJAC. */
7275 /* NDIMEN: Dimension of the space. */
7276 /* MINDGU: Lower limit of the index by U of PATJAC coeff to be taken into account. */
7277 /* MAXDGU: Upper limit of the index by U of PATJAC coeff to be taken into account. */
7278 /* MINDGV: Lower limit of the index by V of PATJAC coeff to be taken into account. */
7279 /* MAXDGV: Upper limit of the index by V of PATJAC coeff to be taken into account. */
7280 /* IORDRU: Order of continuity by U provided by square PATJAC (from -1 to 2) */
7281 /* IORDRV: Order of continuity by V provided by square PATJAC (from -1 to 2) */
7282 /* PATJAC: Table of coeff. of the approximation square with */
7283 /* constraints of order IORDRU by U and IORDRV by V. */
7285 /* OUTPUT ARGUMENTS : */
7286 /* ------------------- */
7287 /* ERRMOY: Average error commited by preserving only the coeff of */
7288 /* PATJAC 2*(IORDRU+1) in MINDGU by U and 2*(IORDRV+1) in MINDGV by V. */
7290 /* COMMONS USED : */
7291 /* ---------------- */
7293 /* REFERENCES CALLED : */
7294 /* --------------------- */
7296 /* DESCRIPTION/NOTES/LIMITATIONS : */
7297 /* ----------------------------------- */
7298 /* Table PATJAC stores the coeff. Cij of */
7299 /* approximation square F(U,V). Indexes i and j show the degree by */
7300 /* U and by V of the base polynoms. These base polynoms are in the form: */
7302 /* ((1 - U*U)**(IORDRU+1)).J(i-2*(IORDRU+1)(U), where */
7304 /* polynom J(i-2*(IORDU+1)(U) is the Jacobi polynom of order */
7305 /* IORDRU+1 (the same by V by replacing U by V in the above expression). */
7307 /* The contribution to the average error of term Cij when */
7308 /* it is removed from PATJAC is Cij*Cij. */
7311 /* ***********************************************************************
7313 /* Name of the routine */
7316 /* ----------------------------- Initialisations ------------------------
7319 /* Parameter adjustments */
7320 patjac_dim1 = *ndgumx + 1;
7321 patjac_dim2 = *ndgvmx + 1;
7322 patjac_offset = patjac_dim1 * patjac_dim2;
7323 patjac -= patjac_offset;
7326 ldbg = AdvApp2Var_SysBase::mnfndeb_() >= 3;
7328 AdvApp2Var_SysBase::mgenmsg_("MMA2MOY", 7L);
7331 idebu = (*iordru + 1) << 1;
7332 idebv = (*iordrv + 1) << 1;
7333 minu = advapp_max(idebu,*mindgu);
7334 minv = advapp_max(idebv,*mindgv);
7338 /* ------------------ Calculation of the upper bound of the average error ------------ */
7339 /* -------------------- when the coeff. of indexes from MINDGU to MAXDGU ------ */
7340 /* ---------------- by U and of indexes from MINDGV to MAXDGV by V are removed -------------- */
7343 for (nd = 1; nd <= i__1; ++nd) {
7345 for (jj = minv; jj <= i__2; ++jj) {
7347 for (ii = idebu; ii <= i__3; ++ii) {
7348 bid1 = patjac[ii + (jj + nd * patjac_dim2) * patjac_dim1];
7349 bid0 += bid1 * bid1;
7358 for (nd = 1; nd <= i__1; ++nd) {
7360 for (jj = idebv; jj <= i__2; ++jj) {
7362 for (ii = minu; ii <= i__3; ++ii) {
7363 bid1 = patjac[ii + (jj + nd * patjac_dim2) * patjac_dim1];
7364 bid0 += bid1 * bid1;
7372 /* ----------------------- Calculation of the average error -------------
7376 *errmoy = sqrt(bid0);
7378 /* ------------------------- The end ------------------------------------
7382 AdvApp2Var_SysBase::mgsomsg_("MMA2MOY", 7L);
7387 //=======================================================================
7388 //function : mma2roo_
7390 //=======================================================================
7391 int AdvApp2Var_ApproxF2var::mma2roo_(integer *nbpntu,
7396 /* System generated locals */
7399 /* Local variables */
7400 static integer ii, ibb;
7402 /* **********************************************************************
7407 /* Return roots of Legendre for discretisations. */
7411 /* TOUS, AB_SPECIFI::CONTRAINTE&, DISCRETISATION, &POINT */
7413 /* INPUT ARGUMENTS : */
7414 /* ------------------ */
7415 /* NBPNTU: Nb of INTERNAL parameters of discretization BY U. */
7416 /* This is also the nb of root of the Legendre polynom where the discretization is done. */
7417 /* NBPNTV: Nb of INTERNAL parameters of discretization BY V. */
7418 /* This is also the nb of root of the Legendre polynom where the discretization is done. */
7420 /* OUTPUT ARGUMENTS : */
7421 /* ------------------- */
7422 /* UROOTL: Table of parameters of discretisation ON (-1,1) BY U.
7424 /* VROOTL: Table of parameters of discretisation ON (-1,1) BY V.
7427 /* COMMONS USED : */
7428 /* ---------------- */
7430 /* REFERENCES CALLED : */
7431 /* --------------------- */
7433 /* DESCRIPTION/NOTES/LIMITATIONS : */
7434 /* ----------------------------------- */
7437 /* **********************************************************************
7440 /* Name of the routine */
7443 /* Parameter adjustments */
7448 ibb = AdvApp2Var_SysBase::mnfndeb_();
7450 AdvApp2Var_SysBase::mgenmsg_("MMA2ROO", 7L);
7453 /* ---------------- Return the POSITIVE roots on U ------------------
7456 AdvApp2Var_MathBase::mmrtptt_(nbpntu, &urootl[(*nbpntu + 1) / 2 + 1]);
7458 for (ii = 1; ii <= i__1; ++ii) {
7459 urootl[ii] = -urootl[*nbpntu - ii + 1];
7462 if (*nbpntu % 2 == 1) {
7463 urootl[*nbpntu / 2 + 1] = 0.;
7466 /* ---------------- Return the POSITIVE roots on V ------------------
7469 AdvApp2Var_MathBase::mmrtptt_(nbpntv, &vrootl[(*nbpntv + 1) / 2 + 1]);
7471 for (ii = 1; ii <= i__1; ++ii) {
7472 vrootl[ii] = -vrootl[*nbpntv - ii + 1];
7475 if (*nbpntv % 2 == 1) {
7476 vrootl[*nbpntv / 2 + 1] = 0.;
7479 /* ------------------------------ The End -------------------------------
7483 AdvApp2Var_SysBase::mgsomsg_("MMA2ROO", 7L);
7487 //=======================================================================
7488 //function : mmmapcoe_
7490 //=======================================================================
7491 int mmmapcoe_(integer *ndim,
7501 /* System generated locals */
7502 integer somtab_dim1, somtab_offset, diftab_dim1, diftab_offset,
7503 crvjac_dim1, crvjac_offset, gsstab_dim1, i__1, i__2, i__3;
7505 /* Local variables */
7506 static integer igss, ikdeb;
7507 static doublereal bidon;
7508 static integer nd, ik, ir, nbroot, ibb;
7510 /* **********************************************************************
7515 /* Calculate the coefficients of polinomial approximation curve */
7516 /* of degree NDGJAC by the method of smallest squares starting from */
7517 /* the discretization of function on the roots of Legendre polynom */
7518 /* of degree NBPNTS. */
7522 /* FONCTION,APPROXIMATION,COEFFICIENT,POLYNOME */
7524 /* INPUT ARGUMENTS : */
7525 /* ------------------ */
7526 /* NDIM : Dimension of the space. */
7527 /* NDGJAC : Max Degree of the polynom of approximation. */
7528 /* The representation in the orthogonal base starts from degree */
7529 /* 0 to degree NDGJAC-2*(JORDRE+1). The polynomial base */
7530 /* is the base of Jacobi of order -1 (Legendre), 0, 1 and 2 */
7531 /* IORDRE : Order of the base of Jacobi (-1,0,1 or 2). Corresponds */
7532 /* to step of constraints, C0,C1 or C2. */
7533 /* NBPNTS : Degree of the polynom of Legendre on the roots which of */
7534 /* are calculated the coefficients of integration by */
7535 /* Gauss method. It is required to set NBPNTS=30,40,50 or 61 */
7536 /* and NDGJAC < NBPNTS. */
7537 /* SOMTAB : Table of F(ti)+F(-ti) with ti in ROOTAB. */
7538 /* DIFTAB : Table of F(ti)-F(-ti) with ti in ROOTAB. */
7539 /* GSSTAB(i,k) : Table of coefficients of integration by the Gauss method : */
7540 /* i varies from 0 to NBPNTS and */
7541 /* k varies from 0 to NDGJAC-2*(JORDRE+1). */
7543 /* OUTPUT ARGUMENTSE : */
7544 /* ------------------- */
7545 /* CRVJAC : Curve of approximation of FONCNP with eventually */
7546 /* taking into account of constraints at the extremities. */
7547 /* This curve is of degree NDGJAC. */
7549 /* COMMONS USED : */
7550 /* ---------------- */
7552 /* REFERENCES CALLED : */
7553 /* --------------------- */
7555 /* DESCRIPTION/NOTES/LIMITATIONS : */
7556 /* ------------------------------- */
7558 /* **********************************************************************
7561 /* Name of the routine */
7563 /* Parameter adjustments */
7564 crvjac_dim1 = *ndgjac + 1;
7565 crvjac_offset = crvjac_dim1;
7566 crvjac -= crvjac_offset;
7567 gsstab_dim1 = *nbpnts / 2 + 1;
7568 diftab_dim1 = *nbpnts / 2 + 1;
7569 diftab_offset = diftab_dim1;
7570 diftab -= diftab_offset;
7571 somtab_dim1 = *nbpnts / 2 + 1;
7572 somtab_offset = somtab_dim1;
7573 somtab -= somtab_offset;
7576 ibb = AdvApp2Var_SysBase::mnfndeb_();
7578 AdvApp2Var_SysBase::mgenmsg_("MMMAPCO", 7L);
7580 ikdeb = (*iordre + 1) << 1;
7581 nbroot = *nbpnts / 2;
7584 for (nd = 1; nd <= i__1; ++nd) {
7586 /* ----------------- Calculate the coefficients of even degree ----------
7590 for (ik = ikdeb; ik <= i__2; ik += 2) {
7594 for (ir = 1; ir <= i__3; ++ir) {
7595 bidon += somtab[ir + nd * somtab_dim1] * gsstab[ir + igss *
7599 crvjac[ik + nd * crvjac_dim1] = bidon;
7603 /* --------------- Calculate the coefficients of uneven degree ----------
7607 for (ik = ikdeb + 1; ik <= i__2; ik += 2) {
7611 for (ir = 1; ir <= i__3; ++ir) {
7612 bidon += diftab[ir + nd * diftab_dim1] * gsstab[ir + igss *
7616 crvjac[ik + nd * crvjac_dim1] = bidon;
7623 /* ------- Add terms connected to the supplementary root (0.D0) ------ */
7624 /* ----------- of Legendre polynom of uneven degree NBPNTS -----------
7627 if (*nbpnts % 2 == 0) {
7631 for (nd = 1; nd <= i__1; ++nd) {
7633 for (ik = ikdeb; ik <= i__2; ik += 2) {
7635 crvjac[ik + nd * crvjac_dim1] += somtab[nd * somtab_dim1] *
7636 gsstab[igss * gsstab_dim1];
7642 /* ------------------------------ The end -------------------------------
7647 AdvApp2Var_SysBase::mgsomsg_("MMMAPCO", 7L);
7651 //=======================================================================
7652 //function : mmaperm_
7654 //=======================================================================
7655 int mmaperm_(integer *ncofmx,
7663 /* System generated locals */
7664 integer crvjac_dim1, crvjac_offset, i__1, i__2;
7666 /* Local variables */
7667 static doublereal bidj;
7668 static integer i__, ia, nd, ncfcut, ibb;
7669 static doublereal bid;
7671 /* **********************************************************************
7676 /* Calculate the square root of the average quadratic error */
7677 /* of approximation done when only the */
7678 /* first NCFNEW coefficients of a curve of degree NCOEFF-1 */
7679 /* written in NORMALIZED Jacobi base of order 2*(IORDRE+1) are preserved. */
7683 /* LEGENDRE,POLYGONE,APPROXIMATION,ERREUR. */
7685 /* INPUT ARGUMENTS : */
7686 /* ------------------ */
7687 /* NCOFMX : Maximum degree of the curve. */
7688 /* NDIM : Dimension of the space. */
7689 /* NCOEFF : Degree +1 of the curve. */
7690 /* IORDRE : Order of constraint of continuity at the extremities. */
7691 /* CRVJAC : The curve the degree which of will be lowered. */
7692 /* NCFNEW : Degree +1 of the resulting polynom. */
7694 /* OUTPUT ARGUMENTS : */
7695 /* ------------------- */
7696 /* ERRMOY : Average precision of approximation. */
7698 /* COMMONS USED : */
7699 /* ---------------- */
7701 /* REFERENCES CALLED : */
7702 /* ----------------------- */
7704 /* DESCRIPTION/NOTES/LIMITATIONS : */
7705 /* ----------------------------------- */
7707 /* ***********************************************************************
7710 /* Name of the routine */
7712 /* Parameter adjustments */
7713 crvjac_dim1 = *ncofmx;
7714 crvjac_offset = crvjac_dim1 + 1;
7715 crvjac -= crvjac_offset;
7718 ibb = AdvApp2Var_SysBase::mnfndeb_();
7720 AdvApp2Var_SysBase::mgenmsg_("MMAPERM", 7L);
7723 /* --------- Minimum degree that can be reached : Stop at 1 or IA -------
7726 ia = (*iordre + 1) << 1;
7728 if (*ncfnew + 1 > ncfcut) {
7729 ncfcut = *ncfnew + 1;
7732 /* -------------- Elimination of coefficients of high degree ------------ */
7733 /* ----------- Loop on the series of Jacobi :NCFCUT --> NCOEFF --------- */
7738 for (nd = 1; nd <= i__1; ++nd) {
7740 for (i__ = ncfcut; i__ <= i__2; ++i__) {
7741 bidj = crvjac[i__ + nd * crvjac_dim1];
7748 /* ----------- Square Root of average quadratic error e -----------
7752 *errmoy = sqrt(bid);
7754 /* ------------------------------- The end ------------------------------
7758 AdvApp2Var_SysBase::mgsomsg_("MMAPERM", 7L);
7762 //=======================================================================
7763 //function : mmapptt_
7765 //=======================================================================
7766 int AdvApp2Var_ApproxF2var::mmapptt_(const integer *ndgjac,
7767 const integer *nbpnts,
7768 const integer *jordre,
7772 /* System generated locals */
7773 integer cgauss_dim1, i__1;
7775 /* Local variables */
7776 static integer kjac, iptt, ipdb0, infdg, iptdb, mxjac, ilong, ibb;
7778 /* **********************************************************************
7783 /* Load the elements required for integration by */
7784 /* Gauss method to obtain the coefficients in the base of
7785 /* Legendre of the approximation by the least squares of a */
7786 /* function. The elements are stored in commons MMAPGSS */
7787 /* (case without constraint), MMAPGS0 (constraints C0), MMAPGS1 */
7788 /* (constraints C1) and MMAPGS2 (constraints C2). */
7792 /* INTEGRATION,GAUSS,JACOBI */
7794 /* INPUT ARGUMENTS : */
7795 /* ------------------ */
7796 /* NDGJAC : Max degree of the polynom of approximation. */
7797 /* The representation in orthogonal base goes from degree
7798 /* 0 to degree NDGJAC-2*(JORDRE+1). The polynomial base */
7799 /* is the base of Jacobi of order -1 (Legendre), 0, 1 and 2 */
7800 /* NBPNTS : Degree of the polynom of Legendre on the roots which of */
7801 /* are calculated the coefficients of integration by the */
7802 /* method of Gauss. It is required that NBPNTS=8,10,15,20,25, */
7803 /* 30,40,50 or 61 and NDGJAC < NBPNTS. */
7804 /* JORDRE : Order of the base of Jacobi (-1,0,1 or 2). Corresponds */
7805 /* to step of constraints C0,C1 or C2. */
7807 /* OUTPUT ARGUMENTS : */
7808 /* ------------------- */
7809 /* CGAUSS(i,k) : Table of coefficients of integration by */
7810 /* Gauss method : i varies from 0 to the integer part */
7811 /* of NBPNTS/2 and k varies from 0 to NDGJAC-2*(JORDRE+1). */
7812 /* These are the coeff. of integration associated to */
7813 /* positive roots of the polynom of Legendre of degree */
7814 /* NBPNTS. CGAUSS(0,k) contains coeff. */
7815 /* of integration associated to root t = 0 when */
7816 /* NBPNTS is uneven. */
7817 /* IERCOD : Error code. */
7819 /* = 11 NBPNTS is not 8,10,15,20,25,30,40,50 or 61. */
7820 /* = 21 JORDRE is not -1,0,1 or 2. */
7821 /* = 31 NDGJAC is too great or too small. */
7823 /* COMMONS USED : */
7824 /* ---------------- */
7825 /* MMAPGSS,MMAPGS0,MMAPGS1,MMAPGS2. */
7826 /* ***********************************************************************
7828 /* Parameter adjustments */
7829 cgauss_dim1 = *nbpnts / 2 + 1;
7832 ibb = AdvApp2Var_SysBase::mnfndeb_();
7834 AdvApp2Var_SysBase::mgenmsg_("MMAPPTT", 7L);
7838 /* ------------------- Tests on the validity of inputs ----------------
7841 infdg = (*jordre + 1) << 1;
7842 if (*nbpnts != 8 && *nbpnts != 10 && *nbpnts != 15 && *nbpnts != 20 && *
7843 nbpnts != 25 && *nbpnts != 30 && *nbpnts != 40 && *nbpnts != 50 &&
7848 if (*jordre < -1 || *jordre > 2) {
7852 if (*ndgjac >= *nbpnts || *ndgjac < infdg) {
7856 /* --------------- Calculation of the start pointer following NBPNTS -----------
7861 iptdb += (8 - infdg) << 2;
7864 iptdb += (10 - infdg) * 5;
7867 iptdb += (15 - infdg) * 7;
7870 iptdb += (20 - infdg) * 10;
7873 iptdb += (25 - infdg) * 12;
7876 iptdb += (30 - infdg) * 15;
7879 iptdb += (40 - infdg) * 20;
7882 iptdb += (50 - infdg) * 25;
7887 ipdb0 = ipdb0 + (14 - infdg) / 2 + 1;
7890 ipdb0 = ipdb0 + (24 - infdg) / 2 + 1;
7893 /* ------------------ Choice of the common depending on JORDRE -------------
7896 if (*jordre == -1) {
7909 /* ---------------- Common MMAPGSS (case without constraints) ----------------
7913 ilong = *nbpnts / 2 << 3;
7915 for (kjac = 0; kjac <= i__1; ++kjac) {
7916 iptt = iptdb + kjac * (*nbpnts / 2) + 1;
7917 AdvApp2Var_SysBase::mcrfill_(&ilong,
7918 &mmapgss_.gslxjs[iptt - 1],
7919 &cgauss[kjac * cgauss_dim1 + 1]);
7922 /* --> Case when the number of points is uneven. */
7923 if (*nbpnts % 2 == 1) {
7926 for (kjac = 0; kjac <= i__1; kjac += 2) {
7927 cgauss[kjac * cgauss_dim1] = mmapgss_.gsl0js[iptt - 1];
7932 for (kjac = 1; kjac <= i__1; kjac += 2) {
7933 cgauss[kjac * cgauss_dim1] = 0.;
7939 /* ---------------- Common MMAPGS0 (case with constraints C0) -------------
7943 mxjac = *ndgjac - infdg;
7944 ilong = *nbpnts / 2 << 3;
7946 for (kjac = 0; kjac <= i__1; ++kjac) {
7947 iptt = iptdb + kjac * (*nbpnts / 2) + 1;
7948 AdvApp2Var_SysBase::mcrfill_(&ilong,
7949 &mmapgs0_.gslxj0[iptt - 1],
7950 &cgauss[kjac * cgauss_dim1 + 1]);
7953 /* --> Case when the number of points is uneven. */
7954 if (*nbpnts % 2 == 1) {
7957 for (kjac = 0; kjac <= i__1; kjac += 2) {
7958 cgauss[kjac * cgauss_dim1] = mmapgs0_.gsl0j0[iptt - 1];
7963 for (kjac = 1; kjac <= i__1; kjac += 2) {
7964 cgauss[kjac * cgauss_dim1] = 0.;
7970 /* ---------------- Common MMAPGS1 (case with constraints C1) -------------
7974 mxjac = *ndgjac - infdg;
7975 ilong = *nbpnts / 2 << 3;
7977 for (kjac = 0; kjac <= i__1; ++kjac) {
7978 iptt = iptdb + kjac * (*nbpnts / 2) + 1;
7979 AdvApp2Var_SysBase::mcrfill_(&ilong,
7980 &mmapgs1_.gslxj1[iptt - 1],
7981 &cgauss[kjac * cgauss_dim1 + 1]);
7984 /* --> Case when the number of points is uneven. */
7985 if (*nbpnts % 2 == 1) {
7988 for (kjac = 0; kjac <= i__1; kjac += 2) {
7989 cgauss[kjac * cgauss_dim1] = mmapgs1_.gsl0j1[iptt - 1];
7994 for (kjac = 1; kjac <= i__1; kjac += 2) {
7995 cgauss[kjac * cgauss_dim1] = 0.;
8001 /* ---------------- Common MMAPGS2 (case with constraints C2) -------------
8005 mxjac = *ndgjac - infdg;
8006 ilong = *nbpnts / 2 << 3;
8008 for (kjac = 0; kjac <= i__1; ++kjac) {
8009 iptt = iptdb + kjac * (*nbpnts / 2) + 1;
8010 AdvApp2Var_SysBase::mcrfill_(&ilong,
8011 &mmapgs2_.gslxj2[iptt - 1],
8012 &cgauss[kjac * cgauss_dim1 + 1]);
8015 /* --> Cas of uneven number of points. */
8016 if (*nbpnts % 2 == 1) {
8019 for (kjac = 0; kjac <= i__1; kjac += 2) {
8020 cgauss[kjac * cgauss_dim1] = mmapgs2_.gsl0j2[iptt - 1];
8025 for (kjac = 1; kjac <= i__1; kjac += 2) {
8026 cgauss[kjac * cgauss_dim1] = 0.;
8032 /* ------------------------- Return the error code --------------
8034 /* --> NBPNTS is not OK */
8038 /* --> JORDRE is not OK */
8042 /* --> NDGJAC is not OK */
8047 /* -------------------------------- The end -----------------------------
8052 AdvApp2Var_SysBase::maermsg_("MMAPPTT", iercod, 7L);
8055 AdvApp2Var_SysBase::mgsomsg_("MMAPPTT", 7L);
8061 //=======================================================================
8062 //function : mmjacpt_
8064 //=======================================================================
8065 int mmjacpt_(const integer *ndimen,
8066 const integer *ncoefu,
8067 const integer *ncoefv,
8068 const integer *iordru,
8069 const integer *iordrv,
8070 const doublereal *ptclgd,
8074 /* System generated locals */
8075 integer ptccan_dim1, ptccan_dim2, ptccan_offset, ptclgd_dim1, ptclgd_dim2,
8076 ptclgd_offset, ptcaux_dim1, ptcaux_dim2, ptcaux_dim3,
8077 ptcaux_offset, i__1, i__2, i__3;
8079 /* Local variables */
8080 static integer kdim, nd, ii, jj, ibb;
8082 /* ***********************************************************************
8087 /* Passage from canonical to Jacobi base for a */
8088 /* "square" in a space of arbitrary dimension. */
8092 /* SMOOTHING,BASE,LEGENDRE */
8095 /* INPUT ARGUMENTS : */
8096 /* ------------------ */
8097 /* NDIMEN : Dimension of the space. */
8098 /* NCOEFU : Degree+1 by U. */
8099 /* NCOEFV : Degree+1 by V. */
8100 /* IORDRU : Order of Jacobi polynoms by U. */
8101 /* IORDRV : Order of Jacobi polynoms by V. */
8102 /* PTCLGD : The square in the Jacobi base. */
8104 /* OUTPUT ARGUMENTS : */
8105 /* ------------------- */
8106 /* PTCAUX : Auxilliary space. */
8107 /* PTCCAN : The square in the canonic base (-1,1) */
8109 /* COMMONS USED : */
8110 /* ---------------- */
8112 /* APPLIED REFERENCES : */
8113 /* ----------------------- */
8115 /* DESCRIPTION/NOTES/LIMITATIONS : */
8116 /* ----------------------------------- */
8117 /* Cancels and replaces MJACPC */
8119 /* *********************************************************************
8121 /* Name of the routine */
8124 /* Parameter adjustments */
8125 ptccan_dim1 = *ncoefu;
8126 ptccan_dim2 = *ncoefv;
8127 ptccan_offset = ptccan_dim1 * (ptccan_dim2 + 1) + 1;
8128 ptccan -= ptccan_offset;
8129 ptcaux_dim1 = *ncoefv;
8130 ptcaux_dim2 = *ncoefu;
8131 ptcaux_dim3 = *ndimen;
8132 ptcaux_offset = ptcaux_dim1 * (ptcaux_dim2 * (ptcaux_dim3 + 1) + 1) + 1;
8133 ptcaux -= ptcaux_offset;
8134 ptclgd_dim1 = *ncoefu;
8135 ptclgd_dim2 = *ncoefv;
8136 ptclgd_offset = ptclgd_dim1 * (ptclgd_dim2 + 1) + 1;
8137 ptclgd -= ptclgd_offset;
8140 ibb = AdvApp2Var_SysBase::mnfndeb_();
8142 AdvApp2Var_SysBase::mgenmsg_("MMJACPT", 7L);
8145 /* Passage into canonical by u. */
8147 kdim = *ndimen * *ncoefv;
8148 AdvApp2Var_MathBase::mmjaccv_(ncoefu,
8151 &ptclgd[ptclgd_offset],
8152 &ptcaux[ptcaux_offset],
8153 &ptccan[ptccan_offset]);
8155 /* Swapping of u and v. */
8158 for (nd = 1; nd <= i__1; ++nd) {
8160 for (jj = 1; jj <= i__2; ++jj) {
8162 for (ii = 1; ii <= i__3; ++ii) {
8163 ptcaux[jj + (ii + (nd + ptcaux_dim3) * ptcaux_dim2) *
8164 ptcaux_dim1] = ptccan[ii + (jj + nd * ptccan_dim2) *
8173 /* Passage into canonical by v. */
8175 kdim = *ndimen * *ncoefu;
8176 AdvApp2Var_MathBase::mmjaccv_(ncoefv,
8179 &ptcaux[((ptcaux_dim3 + 1) * ptcaux_dim2 + 1) * ptcaux_dim1 + 1],
8180 &ptccan[ptccan_offset],
8181 &ptcaux[(((ptcaux_dim3 << 1) + 1) * ptcaux_dim2 + 1) * ptcaux_dim1 + 1]);
8183 /* Swapping of u and v. */
8186 for (nd = 1; nd <= i__1; ++nd) {
8188 for (jj = 1; jj <= i__2; ++jj) {
8190 for (ii = 1; ii <= i__3; ++ii) {
8191 ptccan[ii + (jj + nd * ptccan_dim2) * ptccan_dim1] = ptcaux[
8192 jj + (ii + (nd + (ptcaux_dim3 << 1)) * ptcaux_dim2) *
8201 /* ---------------------------- THAT'S ALL FOLKS ------------------------
8205 AdvApp2Var_SysBase::mgsomsg_("MMJACPT", 7L);