* fichier : pente.dgibi *********************************************************** **** Finite Volume, "Cell-Centred Formulation". **** **** PENT, operator to compute gradients and limiters **** **** **** **** A. BECCANTINI, TTMF JANUARY 2001 **** **** Modif, 10/07/01, syntaxe de PENT changée **** **** Modif, 25/02/04, ajoute 'MODE' 'AXIS' **** *********************************************************** 'ELEM' 'TRI3' 'TRAC' 'X' ; GRAPH = FAUX ; * GRAPH = VRAI ; COMPLET = VRAI ; * *** MESH * P0 = 0.0D0 0.0D0 ; P1 = 3.0D0 0.0D0 ; P2 = 3.0D0 3.0D0 ; P3 = 0.0D0 3.0D0 ; P4 = 6.0D0 0.0D0 ; P5 = 6.0D0 3.0D0 ; 'SI' COMPLET ; N1 = 80 ; N2 = 60 ; N3 = 40 ; N4 = 90 ; N5 = 90 ; N6 = 10 ; N7 = 50 ; 'SINON' ; N1 = 8 ; N2 = 7 ; N3 = 6 ; N4 = 9 ; N5 = 9 ; N6 = 7 ; N7 = 8 ; 'FINSI' ; LINEXT1 = ((P0 'DROIT' N1 P1) 'ET' (P1 'DROIT' N2 P2) 'ET' (P2 'DROIT' N3 P3) 'ET' (P3 'DROIT' N4 P0)) ; LINEXT2 = ((P1 'DROIT' N2 P2) 'ET' (P2 'DROIT' N5 P5) 'ET' (P5 'DROIT' N6 P4) 'ET' (P4 'DROIT' N7 P1)) ; 'OPTION' 'ELEM' QUA4 ; DOM1 = 'SURFACE' LINEXT1 'PLAN' ; 'OPTION' 'ELEM' TRI3 ; DOM2 = 'SURFACE' LINEXT2 'PLAN' ; DOMTOT = DOM1 'ET' DOM2; 'ELIMINATION' 1D-6 DOMTOT; 'MESSAGE' ; 'MESSAGE' 'Le maillage :'; 'MESSAGE' ('CHAINE' 'Nombre TRI3: ' ('EXTRAIRE' LELEM 1)); 'MESSAGE' ('CHAINE' 'Nombre QUA4: ' ('EXTRAIRE' LELEM 2)); * **** Internal and external cells * DOMLIM = 'CONTOUR' DOMTOT 'COULEUR' 'VERTE' ; MDOMTOT = TDOMTOT . 'QUAF' ; MDOMINT = TDOMINT . 'QUAF' ; MDOM1 = TDOM1 . 'QUAF' ; 'SI' GRAPH ; 'TRACER' MAICON 'TITRE' 'Elements sur le contour'; 'TRACER' DOMINT 'TITRE' 'Elements internes' ; 'FINSI' ; * **** The interface points (for limit condition) * X1 Y1 = 'COORDONNEE' POIN0; X0 = X1 ; Y0 = Y1 ; X1 Y1 = 'COORDONNEE' POIN0 ; XFAC = (X0 '+' X1) '/' 2 ; YFAC = (Y0 '+' Y1) '/' 2 ; 'SI' (&BLLIM 'EGA' 1) ; GEOLIM = 'MANUEL' 'POI1' PFAC 'BLEU' ; 'SINON' ; GEOLIM = GEOLIM 'ET' ('MANUEL' 'POI1' PFAC 'BLEU') ; 'FINSI' ; 'FIN' BLLIM ; * * 'ELIMINATION' GEOLIM 1.0D-4 ($DOMTOT . 'FACE') ; * 'SI' GRAPH ; 'TRACER' (DOMTOT 'ET' GEOLIM) 'TITRE' 'DOMTOT et points faces aux bords'; 'FINSI' ; *********************************************** ********* CHP: scalar linear field ********* ********* EULESCAL ********* *********************************************** * * (A11 '*' XX) '+' (A12 '*' YY) '+' A0; * A11 = 12.01517 ; A12 = 13.1421 ; A0 = -3.21 ; B11 = 22.01517 ; B12 = 43.1421 ; B0 = -5.21 ; CHP1 = (A11 '*' XX) '+' (A12 '*' YY) '+' A0; CHP2 = (B11 '*' XX) '+' (B12 '*' YY) '+' B0; CHP = CHP1 'ET' CHP2 ; GRCHP LIMCH GRG = 'PENT' $DOMTOT 'CENTRE' 'EULESCAL' 'NOLIMITE' MOTC CHP ; 'LISTE' MOTC ; 2 'P1' 1.0 'P2' 1.0 ; ERRO = 'MAXIMUM' (LIMCH1 '-' LIMCH) 'ABS' ; 'SI' (ERRO > 1.0D-8) ; 'MESSAGE' ; 'MESSAGE' ; 'ERREUR' 5; 'FINSI' ; * * TEST1 : we check that we have the same values if we use * GRG to compute gradients * GRCHP1 LIMCH = 'PENT' $DOMTOT 'CENTRE' 'EULESCAL' 'LIMITEUR' MOTC CHP 'GRADGEO' GRG ; 'LISTE' MOTC ; ERRO = 'MAXIMUM' (GRCHP '-' GRCHP1) 'ABS' ; 'SI' (ERRO > 1.0D-8) ; 'MESSAGE' ; 'MESSAGE' ; 'ERREUR' 5; 'FINSI' ; * * TEST2 : we check that the computation is exact on the internal domain * 'NATU' 'DISCRET' ; 'NATU' 'DISCRET'; 'NATU' 'DISCRET'; 'NATU' 'DISCRET'; GREXA = GR1DX 'ET' GR1DY 'ET' GR2DX 'ET' GR2DY ; ERRCH = GREXA '-' GRCHPI ; MERRCH = ('PSCAL' ERRCH ERRCH LMOT LMOT) '**' 0.5 ; * *** Error graphics (we change the error CHAMPOINT into a CHAMELEM * 'SI' GRAPH ; 'FINSI' ; ERRO = 'MAXIMUM' ERRCH 'ABS' ; 'SI' (ERRO > 1.0D-8) ; 'MESSAGE' ; 'MESSAGE' ; 'ERREUR' 5; 'FINSI' ; * * TEST3 : we impose the linear field values on the border. * In that case error must be "zero" everywhere. * XLIM YLIM = 'COORDONNEE' GEOLIM ; CHPLI1 = (A11 '*' XLIM) '+' (A12 '*' YLIM) '+' A0; CHPLI2 = (B11 '*' XLIM) '+' (B12 '*' YLIM) '+' B0; CHPLI = CHPLI1 'ET' CHPLI2 ; GRCHPL LIMCH GRGL = 'PENT' $DOMTOT 'CENTRE' 'EULESCAL' 'NOLIMITE' MOTC CHP 'CLIM' CHPLI ; GRCHP1 LIMCH = 'PENT' $DOMTOT 'CENTRE' 'EULESCAL' 'NOLIMITE' MOTC CHP 'CLIM' CHPLI 'GRADGEO' GRGL ; ERRO = 'MAXIMUM' (GRCHPL '-' GRCHP1) 'ABS' ; 'SI' (ERRO > 1.0D-8) ; 'MESSAGE' ; 'MESSAGE' ; 'ERREUR' 5; 'FINSI' ; 'NATU' 'DISCRET' ; 'NATU' 'DISCRET'; 'NATU' 'DISCRET'; 'NATU' 'DISCRET'; GREXA = GR1DX 'ET' GR1DY 'ET' GR2DX 'ET' GR2DY ; ERRCH = GREXA '-' GRCHPL ; MERRCH = ('PSCAL' ERRCH ERRCH LMOT LMOT) '**' 0.5 ; 'SI' GRAPH ; 'TRACER' CHM_MER MOD1 'TITR' ('CHAINE' 'EULESCAL avec C.L. : Erreur'); * CHM_GRA = 'KCHA' $DOMTOT 'CHAM' GRCHPL ; * 'TRACER' CHM_GRA MOD1 'TITR' ('CHAINE' 'Gradient'); 'FINSI' ; ERRO = 'MAXIMUM' ERRCH 'ABS' ; 'SI' (ERRO > 1.0D-8) ; 'MESSAGE' ; 'MESSAGE' ; 'ERREUR' 5; 'FINSI' ; * **** TEST4: boudary condition for a scalar field (EULESCAL) * 'OPTION' 'ELEM' 'QUA4' ; L0 = (P0 'DROIT' P3 'DINI' 0.3 'DFIN' 1 ) ; DOMC0 = L0 'TRANSLATION' 1 (-1.0 0) ; DOMC1 = L0 'TRANSLATION' 5 (5.0 0) 'COULEUR' 'ROUGE' ; DOMC2 = DOMC0 'ET' DOMC1 ; 'ELIMINATION' DOMC2 0.0002 ; 'SI' GRAPH ; 'TRACER' (DOMC1 'ET' DOMC2) 'TITRE' 'Traitement de C.L.' ; 'FINSI' ; MDOMC1 = TDOMC1 . 'QUAF'; MDOMC2 = TDOMC2 . 'QUAF'; * * CHP is symmetric with respect to the line L0 * GRCHP LIMCHP GRG = 'PENT' $DOMC2 'CENTRE' GRCHP1 LIMCHP GRG1 = 'PENT' $DOMC1 'CENTRE' MERRCH = ('PSCAL' ERRCH ERRCH LMOT LMOT) '**' 0.5 ; 'SI' GRAPH ; 'TRACER' CHM_MER MOD1 'FINSI' ; ERRO = 'MAXIMUM' ERRCH 'ABS' ; 'SI' (ERRO > 1.0D-8) ; 'MESSAGE' ; 'MESSAGE' ; 'ERREUR' 5; 'FINSI' ; *************************************************** ********* CHP: vectorial linear field ************ ********* EULEVECT **************** *************************************************** * * (A11 '*' XX) '+' (A12 '*' YY) '+' A0; * A11 = 11.3 ; A12 = 0.23 ; A0 = 1.0 ; B11 = 77.0 ; B12 = 0.6 ; B0 = 9.11 ; CHP1 = (A11 '*' XX) '+' (A12 '*' YY) '+' A0; CHP2 = (B11 '*' XX) '+' (B12 '*' YY) '+' B0; CHP = CHP1 'ET' CHP2 ; GRCHP LIMCHP GRG = 'PENT' $DOMTOT 'CENTRE' 'EULEVECT' 'NOLIMITE' MOTC CHP ; 'LISTE' MOTC ; * * TEST1 : we check that we have the same values if we use * GRG to compute gradients * GRCHP1 LIMCHP = 'PENT' $DOMTOT 'CENTRE' CHP 'GRADGEO' GRG ; ERRO = 'MAXIMUM' (GRCHP '-' GRCHP1) 'ABS' ; 'SI' (ERRO > 5.0D-8) ; 'MESSAGE' ; 'MESSAGE' ; 'ERREUR' 5; 'FINSI' ; * * TEST2 : we check that the computation is exact on the internal domain * 'NATU' 'DISCRET' ; 'NATU' 'DISCRET'; 'NATU' 'DISCRET'; 'NATU' 'DISCRET'; GREXA = GR1DX 'ET' GR1DY 'ET' GR2DX 'ET' GR2DY ; ERRCH = GREXA '-' GRCHPI ; MERRCH = ('PSCAL' ERRCH ERRCH LMOT LMOT) '**' 0.5 ; * *** CREATION DE 'MODEL' POUR GRAPHIQUER LE CHAMELEM *** * 'SI' GRAPH ; 'FINSI' ; ERRO = 'MAXIMUM' ERRCH 'ABS' ; 'SI' (ERRO > 5.0D-8) ; 'MESSAGE' ; 'MESSAGE' ; 'ERREUR' 5; 'FINSI' ; * * TEST3 : we impose the linear field values on the border. * In that case error must be "zero" everywhere. * XLIM YLIM = 'COORDONNEE' GEOLIM ; CHPLI1 = (A11 '*' XLIM) '+' (A12 '*' YLIM) '+' A0; CHPLI2 = (B11 '*' XLIM) '+' (B12 '*' YLIM) '+' B0; CHPLI = CHPLI1 'ET' CHPLI2 ; GRCHPL LIMCHP GRGL = 'PENT' $DOMTOT 'CENTRE' CHP 'CLIM' CHPLI ; GRCHP1 LIMCHP = 'PENT' $DOMTOT 'CENTRE' CHP 'CLIM' CHPLI 'GRADGEO' GRGL ; ERRO = 'MAXIMUM' (GRCHPL '-' GRCHP1) 'ABS' ; 'SI' (ERRO > 1.0D-8) ; 'MESSAGE' ; 'MESSAGE' ; 'ERREUR' 5; 'FINSI' ; 'NATU' 'DISCRET' ; 'NATU' 'DISCRET'; 'NATU' 'DISCRET'; 'NATU' 'DISCRET'; GREXA = GR1DX 'ET' GR1DY 'ET' GR2DX 'ET' GR2DY ; ERRCH = GREXA '-' GRCHPL ; MERRCH = ('PSCAL' ERRCH ERRCH LMOT LMOT) '**' 0.5 ; 'SI' GRAPH ; 'TRACER' CHM_MER MOD1 'TITR' ('CHAINE' 'EULEVECT avec C.L.: erreur'); * CHM_GRA = 'KCHA' $DOMTOT 'CHAM' GRCHPL ; * 'TRACER' CHM_GRA MOD1 'TITR' ('CHAINE' 'Gradient'); 'FINSI' ; ERRO = 'MAXIMUM' ERRCH 'ABS' ; 'SI' (ERRO > 1.0D-8) ; 'MESSAGE' ; 'MESSAGE' ; 'ERREUR' 5; 'FINSI' ; * **** TEST4: boudary condition for a vectorial field (EULESCAL) * 'OPTION' 'ELEM' 'QUA4' ; L0 = (P0 'DROIT' P3 'DINI' 0.3 'DFIN' 1 ) ; DOMC0 = L0 'TRANSLATION' 1 (-1.0 0) ; DOMC1 = L0 'TRANSLATION' 5 (5.0 0) 'COULEUR' 'ROUGE' ; DOMC2 = DOMC0 'ET' DOMC1 ; 'ELIMINATION' DOMC2 0.0002 ; 'SI' GRAPH ; 'TRACER' DOMC2 'TITRE' 'Traitement de C.L.' ; 'FINSI' ; MDOMC1 = TDOMC1 . 'QUAF'; MDOMC2 = TDOMC2 . 'QUAF'; * * L0: symmetry line * GRCHP LIMCHP GRG = 'PENT' $DOMC2 'CENTRE' GRCHP1 LIMCHP GRG1 = 'PENT' $DOMC1 'CENTRE' MERRCH = ('PSCAL' ERRCH ERRCH LMOT LMOT) '**' 0.5 ; 'SI' GRAPH ; 'TRACER' CHM_MER MOD1 * CHM_GRA = 'KCHA' $DOMC1 'CHAM' GRCHP1 ; * 'TRACER' CHM_GRA MOD1 'TITR' ('CHAINE' 'Gradient'); 'FINSI' ; ERRO = 'MAXIMUM' ERRCH 'ABS' ; 'SI' (ERRO > 1.0D-8) ; 'MESSAGE' ; 'MESSAGE' ; 'ERREUR' 5; 'FINSI' ; **************************** **** Diamond method ******** **************************** * * 2D * *** MESH * P0 = 0.0D0 0.0D0 ; P1 = 3.0D0 0.0D0 ; P2 = 3.0D0 3.0D0 ; P3 = 0.0D0 3.0D0 ; P4 = 6.0D0 0.0D0 ; P5 = 6.0D0 3.0D0 ; 'SI' COMPLET ; N1 = 80 ; N2 = 60 ; N3 = 40 ; N4 = 90 ; N5 = 90 ; N6 = 10 ; N7 = 50 ; 'SINON' ; N1 = 8 ; N2 = 6 ; N3 = 4 ; N4 = 9 ; N5 = 9 ; N6 = 3 ; N7 = 5 ; 'FINSI' ; LINEXT1 = ((P0 'DROIT' N1 P1) 'ET' (P1 'DROIT' N2 P2) 'ET' (P2 'DROIT' N3 P3) 'ET' (P3 'DROIT' N4 P0)) ; LINEXT2 = ((P1 'DROIT' N2 P2) 'ET' (P2 'DROIT' N5 P5) 'ET' (P5 'DROIT' N6 P4) 'ET' (P4 'DROIT' N7 P1)) ; 'OPTION' 'ELEM' QUA4 ; DOM1 = 'SURFACE' LINEXT1 'PLAN' ; 'OPTION' 'ELEM' TRI3 ; DOM2 = 'SURFACE' LINEXT2 'PLAN' ; DOMTOT = DOM1 'ET' DOM2; 'ELIMINATION' 1D-6 DOMTOT; 'MESSAGE' ; 'MESSAGE' 'Le maillage :'; 'MESSAGE' ('CHAINE' 'Nombre TRI3: ' ('EXTRAIRE' LELEM 1)); 'MESSAGE' ('CHAINE' 'Nombre QUA4: ' ('EXTRAIRE' LELEM 2)); * **** Internal and external cells * DOMLIM = 'CONTOUR' DOMTOT 'COULEUR' 'VERTE' ; MDOMTOT = TDOMTOT . 'QUAF'; 'SI' GRAPH ; 'TRACER' DOMTOT 'TITRE' 'Domaine' ; 'FINSI' ; A11 = 12.01517 ; A12 = 13.1421 ; A0 = -3.21 ; B11 = 22.01517 ; B12 = 43.1421 ; B0 = -5.21 ; CHP1 = (A11 '*' XX) '+' (A12 '*' YY) '+' A0; CHP2 = (B11 '*' XX) '+' (B12 '*' YY) '+' B0; CHP = CHP1 'ET' CHP2 ; 'DIAMANT' CHP ; * * TEST1 : we check that we have the same values if we use * GRG to compute gradients * 'DIAMANT' CHP 'GRADGEO' GRG ; ERRO = 'MAXIMUM' (GRCHP '-' GRCHP1) 'ABS' ; 'SI' (ERRO > 1.0D-8) ; 'MESSAGE' ; 'MESSAGE' ; 'ERREUR' 5; 'FINSI' ; * * TEST2 : we check that the computation is exact on the internal domain * 'NATU' 'DISCRET' ; 'NATU' 'DISCRET'; 'NATU' 'DISCRET'; 'NATU' 'DISCRET'; GREXA = GR1DX 'ET' GR1DY 'ET' GR2DX 'ET' GR2DY ; ERRCH = GREXA '-' GRCHP1 ; MERRCH = ('PSCAL' ERRCH ERRCH LMOT LMOT) '**' 0.5 ; ERRO = 'MAXIMUM' ERRCH 'ABS' ; 'SI' (ERRO > 1.0D-8) ; 'MESSAGE' ; 'MESSAGE' ; 'ERREUR' 5; 'FINSI' ; * * TEST3 : we impose the linear field values on the border. * In that case error must be "zero" everywhere. * XLIM YLIM = 'COORDONNEE' DOMLIM ; CHPLI1 = (A11 '*' XLIM) '+' (A12 '*' YLIM) '+' A0; CHPLI2 = (B11 '*' XLIM) '+' (B12 '*' YLIM) '+' B0; CHPLI = CHPLI1 'ET' CHPLI2 ; 'DIAMANT' CHP 'CLIM' CHPLI ; 'DIAMANT' CHP 'CLIM' CHPLI 'GRADGEO' GRGL ; ERRO = 'MAXIMUM' (GRCHPL '-' GRCHP1) 'ABS' ; 'SI' (ERRO > 1.0D-8) ; 'MESSAGE' ; 'MESSAGE' ; 'ERREUR' 5; 'FINSI' ; 'NATU' 'DISCRET' ; 'NATU' 'DISCRET'; 'NATU' 'DISCRET'; 'NATU' 'DISCRET'; GREXA = GR1DX 'ET' GR1DY 'ET' GR2DX 'ET' GR2DY ; ERRCH = GREXA '-' GRCHPL ; MERRCH = ('PSCAL' ERRCH ERRCH LMOT LMOT) '**' 0.5 ; ERRO = 'MAXIMUM' ERRCH 'ABS' ; 'SI' (ERRO > 1.0D-8) ; 'MESSAGE' ; 'MESSAGE' ; 'ERREUR' 5; 'FINSI' ; **************************** **** Diamond method ******** **************************** * The BC on flux are taken into account * * 2D * *** MESH * 'OPTION' 'ELEM' QUA4 ; P0 = 0.0D0 0.0D0 ; P1 = 3.0D0 0.0D0 ; P2 = 3.0D0 3.0D0 ; P3 = 0.0D0 3.0D0 ; P4 = 6.0D0 0.0D0 ; P5 = 6.0D0 3.0D0 ; 'SI' COMPLET ; N1 = 80 ; N2 = 60 ; N3 = 40 ; N4 = 90 ; N5 = 90 ; N6 = 10 ; N7 = 50 ; 'SINON' ; N1 = 8 ; N2 = 6 ; N3 = 4 ; N4 = 9 ; N5 = 9 ; N6 = 3 ; N7 = 5 ; 'FINSI' ; LINEXT1 = ((P0 'DROIT' N1 P1) 'ET' (P1 'DROIT' N2 P2) 'ET' (P2 'DROIT' N3 P3) 'ET' (P3 'DROIT' N4 P0)) ; LINEXT2 = ((P1 'DROIT' N2 P2) 'ET' (P2 'DROIT' N5 P5) 'ET' (P5 'DROIT' N6 P4) 'ET' (P4 'DROIT' N7 P1)) ; 'OPTION' 'ELEM' QUA4 ; DOM1 = 'SURFACE' LINEXT1 'PLAN' ; 'OPTION' 'ELEM' TRI3 ; DOM2 = 'SURFACE' LINEXT2 'PLAN' ; DOMTOT = DOM1 'ET' DOM2; 'ELIMINATION' 1D-6 DOMTOT; * BC DOMLIM = 'CONTOUR' DOMTOT 'COULEUR' 'VERTE' ; DOMLIM1 = (DOMLIM 'INTERSECTION' LINEXT1) ; DOMLIM2 = (DOMLIM 'INTERSECTION' LINEXT2) ; * MODELS and MESHES * * QDOMTOT = TDOMTOT . 'QUAF' ; QDOMLIM = TDOMLIM . 'QUAF' ; QDOMLIM1 = TDOMLIM1 . 'QUAF' ; QDOMLIM2 = TDOMLIM2 . 'QUAF' ; 'ELIMINATION' (QDOMTOT 'ET' QDOMLIM 'ET' QDOMLIM1 'ET' QDOMLIM2) 1.0D-6 ; 'SI' GRAPH ; 'FINSI' ; * CHAMPOINT A11 = 12.01517 ; A12 = 13.1421 ; A0 = -3.21 ; B11 = 22.01517 ; B12 = 43.1421 ; B0 = -5.21 ; CHP1 = (A11 '*' XX) '+' (A12 '*' YY) '+' A0; CHP2 = (B11 '*' XX) '+' (B12 '*' YY) '+' B0; CHP = CHP1 'ET' CHP2 ; * * TEST1 : Dirichlet BC everywhere * CHP1 = (A11 '*' XXF) '+' (A12 '*' YYF) '+' A0; CHP2 = (B11 '*' XXF) '+' (B12 '*' YYF) '+' B0; CHPL1 = CHP1 'ET' CHP2 ; CHPL2 = 'COPIER' CACCA ; * 'DIAMAN2' MOTC MOTC1 CHP CHPL1 CHPL2 ; 'LISTE' MOTC ; 'LISTE' MOTC1 ; CHP CHPL1 CHPL2 'GRADGEO' GRG ; ERRO = 'MAXIMUM' (GRCHP '-' GRCHP1) 'ABS' ; 'SI' (ERRO > 1.0D-8) ; 'MESSAGE' ; 'MESSAGE' ; 'ERREUR' 5; 'FINSI' ; 'C1DY' A12 'C2DX' B11 'C2DY' B12 'NATU' 'DISCRET' ; ERRCH = GREXA '-' GRCHP ; ERRO = 'MAXIMUM' ERRCH 'ABS' ; 'SI' (ERRO > 1.0D-8) ; 'MESSAGE' ; 'MESSAGE' ; 'ERREUR' 5; 'FINSI' ; * * TEST2 : von Neumann BC everywhere * 'C1DX' A11 'C1DY' A12 'C2DX' B11 'C2DY' B12 ; CHPL1 = CACCA ; CHP CHPL1 CHPL2 ; ERRO = 'MAXIMUM' (GRCHP '-' GREXA) 'ABS' ; 'SI' (ERRO > 1.0D-8) ; 'MESSAGE' ; 'MESSAGE' ; 'ERREUR' 5; 'FINSI' ; * * TEST3 : Dirichlet BC and von Neumann BC * CHP1 = (A11 '*' XXF) '+' (A12 '*' YYF) '+' A0; CHP2 = (B11 '*' XXF) '+' (B12 '*' YYF) '+' B0; CHPL1 = CHP1 'ET' CHP2 ; 'C1DX' A11 'C1DY' A12 'C2DX' B11 'C2DY' B12 ; CHP CHPL1 CHPL2 ; ERRCH = GREXA '-' GRCHP ; ERRO = 'MAXIMUM' ERRCH 'ABS' ; 'SI' (ERRO > 1.0D-8) ; 'MESSAGE' ; 'MESSAGE' ; 'ERREUR' 5; 'FINSI' ; * ***** Elementary mesh * P0 = 0.0D0 0.0D0 ; P1 = 1.0D0 0.0D0 ; P2 = 1.0D0 1.0D0 ; P3 = 0.0D0 1.0D0 ; LINEXT1 = ((P0 'DROIT' 1 P1) 'ET' (P1 'DROIT' 1 P2) 'ET' (P2 'DROIT' 1 P3) 'ET' (P3 'DROIT' 1 P0)) ; 'OPTION' 'ELEM' 'QUA4' ; DOMTOT = 'SURFACE' LINEXT1 'PLAN' ; * BC DOMLIM = 'CONTOUR' DOMTOT 'COULEUR' 'VERTE' ; * MODELS and MESHES * QUAF QDOMTOT = TDOMTOT . 'QUAF' ; QDOMLIM = TDOMLIM . 'QUAF' ; 'ELIMINATION' (QDOMTOT 'ET' QDOMLIM) 1.0D-6 ; * 'SI' GRAPH ; 'FINSI' ; * CHAMPOINT A11 = 12.01517 ; A12 = 13.1421 ; A0 = -3.21 ; B11 = 22.01517 ; B12 = 43.1421 ; B0 = -5.21 ; CHP1 = (A11 '*' XX) '+' (A12 '*' YY) '+' A0; CHP2 = (B11 '*' XX) '+' (B12 '*' YY) '+' B0; CHP = CHP1 'ET' CHP2 ; * * TEST1 : Dirichlet BC everywhere * CHP1 = (A11 '*' XXF) '+' (A12 '*' YYF) '+' A0; CHP2 = (B11 '*' XXF) '+' (B12 '*' YYF) '+' B0; CHPL1 = CHP1 'ET' CHP2 ; CHPL2 = 'COPIER' CACCA ; * CHP CHPL1 CHPL2 ; CHP CHPL1 CHPL2 'GRADGEO' GRG ; ERRO = 'MAXIMUM' (GRCHP '-' GRCHP1) 'ABS' ; 'SI' (ERRO > 1.0D-8) ; 'MESSAGE' ; 'MESSAGE' ; 'ERREUR' 5; 'FINSI' ; 'C1DY' A12 'C2DX' B11 'C2DY' B12 'NATU' 'DISCRET' ; ERRCH = GREXA '-' GRCHP ; ERRO = 'MAXIMUM' ERRCH 'ABS' ; 'SI' (ERRO > 1.0D-8) ; 'MESSAGE' ; 'MESSAGE' ; 'ERREUR' 5; 'FINSI' ; * * TEST2 : von Neumann BC everywhere * 'C1DX' A11 'C1DY' A12 'C2DX' B11 'C2DY' B12 ; CHPL1 = CACCA ; CHP CHPL1 CHPL2 ; ERRO = 'MAXIMUM' (GRCHP '-' GREXA) 'ABS' ; 'SI' (ERRO > 1.0D-8) ; 'MESSAGE' ; 'MESSAGE' ; 'ERREUR' 5; 'FINSI' ; * * TEST3 : Dirichlet BC and von Neumann BC * XXF YYF = 'COORDONNEE' DOMLIM1 ; CHP1 = (A11 '*' XXF) '+' (A12 '*' YYF) '+' A0; CHP2 = (B11 '*' XXF) '+' (B12 '*' YYF) '+' B0; CHPL1 = CHP1 'ET' CHP2 ; 'C1DX' A11 'C1DY' A12 'C2DX' B11 'C2DY' B12 ; CHP CHPL1 CHPL2 ; ERRCH = GREXA '-' GRCHP ; ERRO = 'MAXIMUM' ERRCH 'ABS' ; 'SI' (ERRO > 1.0D-8) ; 'MESSAGE' ; 'MESSAGE' ; 'ERREUR' 5; 'FINSI' ; ************************************ ************************************ ************************************ ******** EULESCAL AXI ************** ************************************ ************************************ ************************************ ************************************ 'ELEM' 'TRI3' 'TRAC' 'X' 'MODE' 'AXIS' ; GRAPH = FAUX ; * GRAPH = VRAI ; COMPLET = VRAI ; * *** MESH * P0 = 1.0D0 0.0D0 ; P1 = 3.0D0 0.0D0 ; P2 = 3.0D0 3.0D0 ; P3 = 1.0D0 3.0D0 ; P4 = 6.0D0 0.0D0 ; P5 = 6.0D0 3.0D0 ; 'SI' COMPLET ; N1 = 80 ; N2 = 60 ; N3 = 40 ; N4 = 90 ; N5 = 90 ; N6 = 10 ; N7 = 50 ; 'SINON' ; N1 = 8 ; N2 = 7 ; N3 = 6 ; N4 = 9 ; N5 = 9 ; N6 = 7 ; N7 = 8 ; 'FINSI' ; LINEXT1 = ((P0 'DROIT' N1 P1) 'ET' (P1 'DROIT' N2 P2) 'ET' (P2 'DROIT' N3 P3) 'ET' (P3 'DROIT' N4 P0)) ; LINEXT2 = ((P1 'DROIT' N2 P2) 'ET' (P2 'DROIT' N5 P5) 'ET' (P5 'DROIT' N6 P4) 'ET' (P4 'DROIT' N7 P1)) ; 'OPTION' 'ELEM' QUA4 ; DOM1 = 'SURFACE' LINEXT1 'PLAN' ; 'OPTION' 'ELEM' TRI3 ; DOM2 = 'SURFACE' LINEXT2 'PLAN' ; DOMTOT = DOM1 'ET' DOM2; 'ELIMINATION' 1D-6 DOMTOT; 'MESSAGE' ; 'MESSAGE' 'Le maillage :'; 'MESSAGE' ('CHAINE' 'Nombre TRI3: ' ('EXTRAIRE' LELEM 1)); 'MESSAGE' ('CHAINE' 'Nombre QUA4: ' ('EXTRAIRE' LELEM 2)); * **** Internal and external cells * DOMLIM = 'CONTOUR' DOMTOT 'COULEUR' 'VERTE' ; MDOMTOT = TDOMTOT . 'QUAF' ; MDOMINT = TDOMINT . 'QUAF' ; MDOM1 = TDOM1 . 'QUAF' ; 'SI' GRAPH ; 'TRACER' MAICON 'TITRE' 'Elements sur le contour'; 'TRACER' DOMINT 'TITRE' 'Elements internes' ; 'FINSI' ; * **** The interface points (for limit condition) * X1 Y1 = 'COORDONNEE' POIN0; X0 = X1 ; Y0 = Y1 ; X1 Y1 = 'COORDONNEE' POIN0 ; XFAC = (X0 '+' X1) '/' 2 ; YFAC = (Y0 '+' Y1) '/' 2 ; 'SI' (&BLLIM 'EGA' 1) ; GEOLIM = 'MANUEL' 'POI1' PFAC 'BLEU' ; 'SINON' ; GEOLIM = GEOLIM 'ET' ('MANUEL' 'POI1' PFAC 'BLEU') ; 'FINSI' ; 'FIN' BLLIM ; * * 'ELIMINATION' GEOLIM 1.0D-4 ($DOMTOT . 'FACE') ; * 'SI' GRAPH ; 'TRACER' (DOMTOT 'ET' GEOLIM) 'TITRE' 'DOMTOT et points faces aux bords'; 'FINSI' ; *********************************************** ********* CHP: scalar linear field ********* ********* EULESCAL ********* *********************************************** * * (A11 '*' XX) '+' (A12 '*' YY) '+' A0; * A11 = 12.01517 ; A12 = 13.1421 ; A0 = -3.21 ; B11 = 22.01517 ; B12 = 43.1421 ; B0 = -5.21 ; CHP1 = (A11 '*' XX) '+' (A12 '*' YY) '+' A0; CHP2 = (B11 '*' XX) '+' (B12 '*' YY) '+' B0; CHP = CHP1 'ET' CHP2 ; * * TEST3 : we impose the linear field values on the border. * In that case error must be "zero" everywhere. * XLIM YLIM = 'COORDONNEE' GEOLIM ; CHPLI1 = (A11 '*' XLIM) '+' (A12 '*' YLIM) '+' A0; CHPLI2 = (B11 '*' XLIM) '+' (B12 '*' YLIM) '+' B0; CHPLI = CHPLI1 'ET' CHPLI2 ; GRCHPL LIMCH GRGL = 'PENT' $DOMTOT 'CENTRE' GRCHP1 LIMCH = 'PENT' $DOMTOT 'CENTRE' CHP 'CLIM' CHPLI 'GRADGEO' GRGL ; ERRO = 'MAXIMUM' (GRCHPL '-' GRCHP1) 'ABS' ; 'SI' (ERRO > 1.0D-8) ; 'MESSAGE' ; 'MESSAGE' ; 'ERREUR' 5; 'FINSI' ; 'NATU' 'DISCRET' ; 'NATU' 'DISCRET'; 'NATU' 'DISCRET'; 'NATU' 'DISCRET'; GREXA = GR1DX 'ET' GR1DY 'ET' GR2DX 'ET' GR2DY ; ERRCH = GREXA '-' GRCHPL ; MERRCH = ('PSCAL' ERRCH ERRCH LMOT LMOT) '**' 0.5 ; 'SI' GRAPH ; 'TRACER' CHM_MER MOD1 'TITR' ('CHAINE' 'EULESCAL avec C.L. : Erreur'); * CHM_GRA = 'KCHA' $DOMTOT 'CHAM' GRCHPL ; * 'TRACER' CHM_GRA MOD1 'TITR' ('CHAINE' 'Gradient'); 'FINSI' ; ERRO = 'MAXIMUM' ERRCH 'ABS' ; 'SI' (ERRO > 1.0D-8) ; 'MESSAGE' ; 'MESSAGE' ; 'ERREUR' 5; 'FINSI' ; ************************************ ************************************ **** Diamond method with BC ******** ************************************ ************************************ * * 3D * *** MESH * * 'MESSAGE' ; 'MESSAGE' '3D' ; 'MESSAGE' 'PENTE DIAMAN2' ; 'MESSAGE' ; P1 = 0.0 0.0 0.0 ; P2 = 3.0 0.0 0.0 ; P3 = 3.0 3.0 0.0 ; P4 = 0.0 3.0 0.0 ; BAS = 'SURFACE' (P1 'DROIT' 4 P2 'DROIT' 3 P3 'DROIT' 5 P4 'DROIT' 6 P1) 'PLAN' ; DOMTOT = BAS 'VOLUME' 'TRANSLATION' 2 (0.0 0.0 2.0) ; CONDOM = 'ENVELOPPE' DOMTOT ; MDOMTOT = 'MODELISER' DOMTOT 'EULER' ; MCONDOM = 'MODELISER' CONDOM 'EULER' ; MBAS = 'MODELISER' BAS 'EULER' ; MRESTEN = 'MODELISER' RESTEN 'EULER' ; QDOMTOT = TDOMTOT . 'QUAF' ; QCONDOM = TCONDOM . 'QUAF' ; QBAS = TBAS . 'QUAF' ; QRESTEN = TRESTEN . 'QUAF' ; 'ELIMINATION' (QDOMTOT 'ET' QCONDOM 'ET' QBAS 'ET' QRESTEN) 1.0D-6 ; 'SI' GRAPH ; 'FINSI' ; * CHAMPOINT A11 = 12.01517 ; A12 = 13.1421 ; A13 = 11.7 ; A0 = -3.21 ; B11 = 22.01517 ; B12 = 43.1421 ; B13 = 7.01 ; B0 = -5.21 ; CHP1 = (A11 '*' XX) '+' (A12 '*' YY) '+' (A13 '*' ZZ) '+' A0; CHP2 = (B11 '*' XX) '+' (B12 '*' YY) '+' (B13 '*' ZZ) '+' B0; CHP = CHP1 'ET' CHP2 ; * * TEST1 : Dirichlet BC everywhere * CHP1 = (A11 '*' XXF) '+' (A12 '*' YYF) '+' (A13 '*' ZZF) '+' A0; CHP2 = (B11 '*' XXF) '+' (B12 '*' YYF) '+' (B13 '*' ZZF) '+' B0; CHPL1 = CHP1 'ET' CHP2 ; CHPL2 = 'COPIER' CACCA ; * 'C2DX' 'C2DY' 'C2DZ') CHP CHPL1 CHPL2 ; 'C2DX' 'C2DY' 'C2DZ') CHP CHPL1 CHPL2 'GRADGEO' GRG ; ERRO = 'MAXIMUM' (GRCHP '-' GRCHP1) 'ABS' ; 'SI' (ERRO > 1.0D-8) ; 'MESSAGE' ; 'MESSAGE' ; 'ERREUR' 5; 'FINSI' ; 'C1DY' A12 'C1DZ' A13 'C2DX' B11 'C2DY' B12 'C2DZ' B13 'NATU' 'DISCRET' ; ERRCH = GREXA '-' GRCHP ; ERRO = 'MAXIMUM' ERRCH 'ABS' ; 'SI' (ERRO > 1.0D-8) ; 'MESSAGE' ; 'MESSAGE' ; 'ERREUR' 5; 'FINSI' ; * * TEST2 : von Neumann BC everywhere * 'C1DX' A11 'C1DY' A12 'C1DZ' A13 'C2DX' B11 'C2DY' B12 'C2DZ' B13 'NATU' 'DISCRET' ; CHPL1 = CACCA ; 'C2DX' 'C2DY' 'C2DZ') CHP CHPL1 CHPL2 ; ERRO = 'MAXIMUM' (GRCHP '-' GREXA) 'ABS' ; 'SI' (ERRO > 1.0D-8) ; 'MESSAGE' ; 'MESSAGE' ; 'ERREUR' 5; 'FINSI' ; * * TEST2 : von Neumann and Dirichlet BC * CHP1 = (A11 '*' XXF) '+' (A12 '*' YYF) '+' (A13 '*' ZZF) '+' A0; CHP2 = (B11 '*' XXF) '+' (B12 '*' YYF) '+' (B13 '*' ZZF) '+' B0; CHPL1 = CHP1 'ET' CHP2 ; 'C1DX' A11 'C1DY' A12 'C1DZ' A13 'C2DX' B11 'C2DY' B12 'C2DZ' B13 'NATU' 'DISCRET' ; 'C2DX' 'C2DY' 'C2DZ') CHP CHPL1 CHPL2 ; ERRO = 'MAXIMUM' (GRCHP '-' GREXA) 'ABS' ; 'SI' (ERRO > 1.0D-8) ; 'MESSAGE' ; 'MESSAGE' ; 'ERREUR' 5; 'FINSI' ; * * ***************************** ***************************** ***** Maillage regulier ***** ***************************** ***************************** * * Géométrie * pA = 0. 0. 0. ; pB = 1. 0. 0. ; pC = 1. 1. 0. ; pD = 0. 1. 0. ; pE = 0. 0. 1. ; pF = 1. 0. 1. ; pG = 1. 1. 1. ; pH = 0. 1. 1. ; * * Paramètres de la discrétisation de base * nAB = 1 ; nBC = 1 ; nCD = 1 ; nDA = 1 ; nH = 1 ; * * Géométrie discrétisée * bas = 'DROIT' nAB pA pB ; droite = 'DROIT' nBC pB pC ; haut = 'DROIT' nCD pC pD ; gauche = 'DROIT' nDA pD pA ; spourt = bas 'ET' droite 'ET' haut 'ET' gauche ; sbas = 'SURFACE' spourt 'PLAN' ; * 'SI' (graph ET FAUX) ; 'TRACER' (mt 'ET' cnt) 'TITRE' 'Enveloppe' ; 'FINSI' ; * * Creation des modèles * * MODELS * QUAF Mmt = tmt . 'QUAF' ; Mcnt = tcnt . 'QUAF' ; * MAILLAGES 'SI' graph ; 'TITRE' 'Centres et faces' ; 'FINSI' ; * * solex = (alpha '*' x) '+' (beta '*' y) '+' (gamma '*' z) '+' * (eta '*' x '*' y '*' z) '+' delta * alpha = (** 2. 0.5) ; beta = (** 3. 0.5) ; gamma = (** 5. 0.5) ; * eta = (** 7. 0.5) ; eta = 0. ; delta = (** 9. 0.5) ; * * Solution exacte bilinéaire (sur le centre du contour) * solex = (xx '*' alpha) '+' (yy '*' beta) '+' (zz '*' gamma) '+' (xx '*' yy '*' zz '*' eta) '+' delta ; * * Solution exacte aux centres * solexc = (xxc '*' alpha) '+' (yyc '*' beta) '+' (zzc '*' gamma) '+' (xxc '*' yyc '*' zzc '*' eta) '+' delta ; * * Gradient exact (aux faces) * 'DISCRET' ; 'DISCRET' ; 'DISCRET' ; gradex = gradtx 'ET' gradty 'ET' gradtz ; * * BC * * * TEST1 : Dirichlet BC everywhere * * cvn = set of FACE centers where we impose von Neumann boundary * conditions 'SI' graph ; 'TRACER' (cnt et cvn) 'TITRE' 'Von Neumann B.C.' ; 'TRACER' (cnt et cdir) 'TITRE' 'Dirichlet B.C.' ; 'FINSI' ; gradtl ; ERRO = 'MAXIMUM' (gradt '-' gradex) 'ABS' ; 'SI' (ERRO > 1.0D-8) ; 'MESSAGE' ; 'MESSAGE' ; 'ERREUR' 5; 'FINSI' ; * * TEST2 : von Neumann BC everywhere * * cvn = set of FACE centers where we impose von Neumann boundary * conditions 'SI' graph ; 'TRACER' (cnt et cvn) 'TITRE' 'Von Neumann B.C.' ; 'TRACER' (cnt et cdir) 'TITRE' 'Dirichlet B.C.' ; 'FINSI' ; gradtl ; ERRO = 'MAXIMUM' (gradt '-' gradex) 'ABS' ; 'SI' (ERRO > 1.0D-8) ; 'MESSAGE' ; 'MESSAGE' ; 'ERREUR' 5; 'FINSI' ; * * TEST3 : von Neumann and Dirichlet bc * * cvn = set of FACE centers where we impose von Neumann boundary * conditions 'SI' graph ; 'TRACER' (cnt et cvn) 'TITRE' 'Von Neumann B.C.' ; 'TRACER' (cnt et cdir) 'TITRE' 'Dirichlet B.C.' ; 'FINSI' ; gradtl ; ERRO = 'MAXIMUM' (gradt '-' gradex) 'ABS' ; 'SI' (ERRO > 1.0D-8) ; 'MESSAGE' ; 'MESSAGE' ; 'ERREUR' 5; 'FINSI' ; 'FIN' ;
© Cast3M 2003 - Tous droits réservés.
Mentions légales