[SCM] Gerris Flow Solver branch, upstream, updated. b3aa46814a06c9cb2912790b23916ffb44f1f203
Stephane Popinet
popinet at users.sf.net
Fri May 15 02:55:46 UTC 2009
The following commit has been merged in the upstream branch:
commit 9872734f7b578230c65096d98bfe95054bfb542f
Author: Stephane Popinet <popinet at users.sf.net>
Date: Thu Oct 2 11:55:47 2008 +1000
New version of the terrain module using region_sum queries
This patch contains both the old and the new version, for debugging purposes.
darcs-hash:20081002015547-d4795-32f84830aea0ced8e81575e26ff0ecfa1b475cbb.gz
diff --git a/modules/cform.lisp b/modules/cform.lisp
new file mode 100644
index 0000000..79f8abf
--- /dev/null
+++ b/modules/cform.lisp
@@ -0,0 +1,198 @@
+;;;
+;;; cform.lisp -- Maxima output formatter for Programming Language C.
+;;; Copyright (C) 2007-2008 Tomohide Naniwa
+;;; version 1.2: Aug. 8, 2008
+;;; based on precious contribution by D.C. Hauagge
+
+;;; cform.lisp is free software; you can redistribute it
+;;; and/or modify it under the terms of the GNU General Public
+;;; License as published by the Free Software Foundation; either
+;;; version 2, or (at your option) any later version.
+
+;;; cform.lisp is distributed in the hope that it will be
+;;; useful, but WITHOUT ANY WARRANTY; without even the implied
+;;; warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
+;;; See the GNU General Public License for more details.
+
+;;; Based on f90.lisp. Copyright statements for f90.lisp follow:
+;;; Copyright (C) 2004 James F. Amundson
+
+;;; Based on fortra.lisp. Copyright statements for fortra.lisp follow:
+;;; Copyright (c) 1984,1987 by William Schelter,University of Texas
+;;; All rights reserved
+;;; (c) Copyright 1980 Massachusetts Institute of Technology
+
+(in-package "MAXIMA")
+(macsyma-module cform)
+
+(DECLARE-TOP (SPECIAL LB RB ;Used for communication with MSTRING.
+ $LOADPRINT ;If NIL, no load message gets printed.
+ 1//2 -1//2)
+ (*LEXPR C-PRINT $CFORMMX))
+
+(DEFMSPEC $CFORM (L)
+ (SETQ L (FEXPRCHECK L))
+ (LET ((VALUE (STRMEVAL L)))
+ (COND ((MSETQP L) (SETQ VALUE `((MEQUAL) ,(CADR L) ,(MEVAL L)))))
+ (COND ((AND (SYMBOLP L) ($MATRIXP VALUE))
+ ($CFORMMX L VALUE))
+ ((AND (NOT (ATOM VALUE)) (EQ (CAAR VALUE) 'MEQUAL)
+ (SYMBOLP (CADR VALUE)) ($MATRIXP (CADDR VALUE)))
+ ($CFORMMX (CADR VALUE) (CADDR VALUE)))
+ (T (C-PRINT VALUE)))))
+
+;; Some aliases for C language may be omittable for Maxima 5.10.x
+(setq c-alias
+ '(($POW "pow")
+ ($EXP "exp")
+ (%SQRT "sqrt")
+ (%SIN "sin")
+ (%COS "cos")
+ (%TAN "tan")
+ (%ACOS "acos")
+ (%ASIN "asin")
+ (%ATAN "atan")
+ ($ATAN2 "atan2")
+ (%LOG "log")
+ (%SINH "sinh")
+ (%COSH "cosh")
+ (%TANH "tanh")
+ (%ASINH "asinh")
+ (%ACOSH "acosh")
+ (%ATANH "atanh")
+ (MABS "fabs")
+ ))
+
+(DEFUN C-PRINT (X &OPTIONAL (STREAM #+Maclisp NIL #-Maclisp *standard-output*)
+ &AUX #+PDP10 (TERPRI T) #+PDP10 ($LOADPRINT NIL)
+ ;; This is a poor way of saying that array references
+ ;; are to be printed with parens instead of brackets.
+ (LB #\[ ) (RB #\] )
+ ;; Definition of heading white space.
+ (WHS " "))
+ ;; Restructure the expression for displaying.
+ (SETQ X (CSCAN X))
+ ;; Linearize the expression using MSTRING. Some global state must be
+ ;; modified for MSTRING to generate using C syntax. This must be
+ ;; undone so as not to modifiy the toplevel behavior of MSTRING.
+ (UNWIND-PROTECT
+ (DEFPROP MEXPT MSIZE-INFIX GRIND)
+ (DEFPROP MMINUS 100. LBP)
+ (DEFPROP MSETQ (#\:) STRSYM)
+ (mapc (lambda (x)
+ (putprop (car x) (get (car x) 'REVERSEALIAS) 'KEEP-RA)
+ (putprop (car x) (cadr x) 'REVERSEALIAS)) c-alias)
+ (SETQ X (mstring x))
+ ;; Make sure this gets done before exiting this frame.
+ (DEFPROP MEXPT MSZ-MEXPT GRIND)
+ (REMPROP 'MMINUS 'LBP)
+ (mapc (lambda (x)
+ (putprop (car x) (get (car x) 'KEEP-RA) 'REVERSEALIAS)
+ (remprop (car x) 'KEEP-RA)) c-alias)
+ )
+
+ ;; MSTRING returns a list of characters. Now print them
+ (do ((char 0 (1+ char))
+ (line ""))
+ ((>= char (length x)))
+ (setf line (concatenate 'string line (make-sequence
+ 'string 1
+ :initial-element (nth char x))))
+ (if (>= (length line) 80)
+ (let ((break_point -1))
+ (mapc #'(lambda (x)
+ (let ((p (search x line :from-end t)))
+ (if (and p (> p 0))
+ (setf break_point p))))
+ '("+" "-" "*" "/"))
+; (increment break_point)
+ (setf break_point (+ break_point 1))
+ (if (= break_point 0)
+ (progn (princ line stream) (setf line WHS))
+ (progn
+ (princ (subseq line 0 break_point) stream)
+ (terpri stream)
+ (setf line (concatenate 'string WHS
+ (subseq line break_point
+ (length line))))))))
+ (if (and (= char (1- (length x))) (not (equal line WHS)))
+ (princ line stream))
+ )
+ (princ ";" stream)
+ (terpri stream)
+ '$done)
+
+(DEFUN CSCAN (E)
+ (COND ((ATOM E) (cond ((eq e '$%i) '((mprogn) 0.0 1.0))
+ (t E))) ;%I is (0,1)
+; Recent C compilers may have prototype declarathions for math functions.
+; ((AND (EQ (CAAR E) 'MEXPT) (EQ (CADR E) '$%E) (numberp (caddr e)))
+; (LIST '($EXP SIMP) (float (CADDR E))))
+ ((AND (EQ (CAAR E) 'MEXPT) (EQ (CADR E) '$%E))
+ (LIST '($EXP SIMP) (CSCAN (CADDR E))))
+ ((AND (EQ (CAAR E) 'MEXPT) (ALIKE1 (CADDR E) 1//2))
+ (LIST '(%SQRT SIMP) (CSCAN (CADR E))))
+ ((AND (EQ (CAAR E) 'MEXPT) (ALIKE1 (CADDR E) -1//2))
+ (LIST '(MQUOTIENT SIMP) 1 (LIST '(%SQRT SIMP) (CSCAN (CADR E)))))
+; ((and (EQ (CAAR E) 'MEXPT) (numberp (caddr E)))
+; (LIST '($POW SIMP) (CSCAN (CADR E)) (float (CADDR E))))
+; ((and (EQ (CAAR E) 'MEXPT) (numberp (cadr E)))
+; (LIST '($POW SIMP) (float (CADR E)) (CSCAN (CADDR E))))
+ ((EQ (CAAR E) 'MEXPT)
+ (LIST '($POW SIMP) (CSCAN (CADR E)) (CSCAN (CADDR E))))
+ ((AND (EQ (CAAR E) 'MTIMES) (RATNUMP (CADR E))
+ (ZL-MEMBER (CADADR E) '(1 -1)))
+ (COND ((EQUAL (CADADR E) 1) (CSCAN-MTIMES E))
+ (T (LIST '(MMINUS SIMP) (CSCAN-MTIMES E)))))
+ ((EQ (CAAR E) 'RAT)
+ (LIST '(MQUOTIENT SIMP) (FLOAT (CADR E)) (FLOAT (CADDR E))))
+ ((EQ (CAAR E) 'MRAT) (CSCAN (RATDISREP E)))
+
+ ;; x[1,2,3] => x[1][2][3]
+ ((AND (EQ (CADDAR E) 'ARRAY) (NOT (EQ (CAAR E) 'MQAPPLY)) )
+ (if (> (LENGTH E) 2)
+ ;; then
+ (LIST '(MQAPPLY SIMP ARRAY)
+ (CSCAN (APPEND (LIST (LIST (CAAR E) 'SIMP 'ARRAY))
+ (BUTLAST (CDR E))))
+ (CAR (LAST E)))
+ ;; else
+ E
+ )
+ )
+
+ ;; complex numbers to f77 syntax a+b%i ==> (a,b)
+ ((and (memq (caar e) '(mtimes mplus))
+ ((lambda (a)
+ (and (numberp (cadr a))
+ (numberp (caddr a))
+ (not (zerop1 (cadr a)))
+ (list '(mprogn) (caddr a) (cadr a))))
+ (simplify ($bothcoef e '$%i)))))
+ (T (CONS (CAR E) (MAPCAR 'CSCAN (CDR E))))))
+
+(DEFUN CSCAN-MTIMES (E)
+ (LIST '(MQUOTIENT SIMP)
+ (COND ((NULL (CDDDR E)) (CSCAN (CADDR E)))
+ (T (CONS (CAR E) (MAPCAR 'CSCAN (CDDR E)))))
+ (FLOAT (CADDR (CADR E)))))
+
+;; Takes a name and a matrix and prints a sequence of C assignment
+;; statements of the form
+;; NAME[I][J] = <corresponding matrix element>
+;; The indcies I, J will be counted from 1.
+
+(DEFMFUN $CFORMMX (NAME MAT &OPTIONAL (STREAM #-CL NIL #+CL *standard-output*)
+ &AUX ($LOADPRINT NIL) (K 'array))
+ (COND ((NOT (symbolp NAME))
+ (MERROR "~%First argument to CFORMMX must be a symbol."))
+ ((NOT ($MATRIXP MAT))
+ (MERROR "Second argument to CFORMMX not a matrix: ~M" MAT)))
+ (DO ((MAT (CDR MAT) (CDR MAT)) (I 1 (f1+ I))) ((NULL MAT))
+ (DECLARE (FIXNUM I))
+ (DO ((M (CDAR MAT) (CDR M)) (J 1 (f1+ J))) ((NULL M))
+ (DECLARE (FIXNUM J))
+ (C-PRINT `((MEQUAL) ((((,NAME ,K) ,I) ,K) ,J) ,(CAR M)) STREAM)))
+ '$DONE)
+
+;; End:
diff --git a/modules/terrain.mac b/modules/terrain.mac
new file mode 100644
index 0000000..ee69b1c
--- /dev/null
+++ b/modules/terrain.mac
@@ -0,0 +1,145 @@
+/* Maxima code to compute coefficients used in terrain.mod:update_terrain() */
+
+/* usage:
+ % maxima -p cform.lisp
+ (%i1) batchload("terrain.mac");
+*/
+
+/* load("cform"); */
+
+hfactor(expr):=
+ block (
+ [a: expand(expr)],
+ a: ratsubst(s_H4,z*z,a),
+ a: ratsubst(s_H3,x*y*z,a),
+ a: expand(a),
+ a: ratsubst(s_m33,x^2*y^2,a),
+ a: expand(a),
+ a: ratsubst(s_m13,x^2*y,a),
+ a: expand(a),
+ a: ratsubst(s_m23,y^2*x,a),
+ a: expand(a),
+ a: ratsubst(s_m11,x^2,a),
+ a: expand(a),
+ a: ratsubst(s_m22,y^2,a),
+ a: expand(a),
+ a: ratsubst(s_m03,x*y,a),
+ a: expand(a),
+ a: ratsubst(s_H1,x*z,a),
+ a: expand(a),
+ a: ratsubst(s_H2,y*z,a),
+ a: expand(a),
+ a: ratsubst(s_m01,x,a),
+ a: expand(a),
+ a: ratsubst(s_m02,y,a),
+ a: expand(a),
+ a: ratsubst(s_H0,z,a),
+ a: expand(a),
+ factorfacsum(a,s_H0,s_H1,s_H2,s_H3,s_H4,s_m01,s_m02,s_m03,s_m13,s_m23,s_m11,s_m22,s_m33)
+ );
+
+/* We want to "recenter" the H coefficient on (xc,yc,hp(x)) where
+ hp(x) is the height on the coarser mesh (i.e. the parent cell) */
+
+H1: (x-xc)*(z-hp[0]-hp[1]*x-hp[2]*y-hp[3]*x*y);
+print("H1=",H1);
+H1: hfactor(H1);
+print("C code:");
+cform(H1);
+print("----------------------------------------------------------------");
+
+H2: (y-yc)*(z-hp[0]-hp[1]*x-hp[2]*y-hp[3]*x*y);
+print("H2=",H2);
+H2: hfactor(H2);
+print("C code:");
+cform(H2);
+print("----------------------------------------------------------------");
+
+H3: (x-xc)*(y-yc)*(z-hp[0]-hp[1]*x-hp[2]*y-hp[3]*x*y);
+print("H3=",H3);
+H3: hfactor(H3);
+print("C code:");
+cform(H3);
+print("----------------------------------------------------------------");
+
+H4: (z-hp[0]-hp[1]*x-hp[2]*y-hp[3]*x*y)^2;
+print("H4=",H4);
+H4: hfactor(H4);
+print("C code:");
+cform(H4);
+print("----------------------------------------------------------------");
+
+/* We want to "recenter" the m coefficient on (xc,yc) */
+
+m03: (x - xc)*(y - yc);
+print("m03=",m03);
+m03: hfactor(m03);
+print("C code:");
+cform(m03);
+print("----------------------------------------------------------------");
+
+m11: (x - xc)*(x - xc);
+print("m11=",m11);
+m11: hfactor(m11);
+print("C code:");
+cform(m11);
+print("----------------------------------------------------------------");
+
+m22: (y - yc)*(y - yc);
+print("m22=",m22);
+m22: hfactor(m22);
+print("C code:");
+cform(m22);
+print("----------------------------------------------------------------");
+
+m13: (x - xc)*(x - xc)*(y - yc);
+print("m13=",m13);
+m13: hfactor(m13);
+print("C code:");
+cform(m13);
+print("----------------------------------------------------------------");
+
+m23: (x - xc)*(y - yc)*(y - yc);
+print("m23=",m23);
+m23: hfactor(m23);
+print("C code:");
+cform(m23);
+print("----------------------------------------------------------------");
+
+m33: (x - xc)*(x - xc)*(y - yc)*(y - yc);
+print("m33=",m33);
+m33: hfactor(m33);
+print("C code:");
+cform(m33);
+print("----------------------------------------------------------------");
+
+/* Only recenter the H coefficients using (xc,yc) if we are looking at
+ the absolute height */
+
+H1: (x-xc)*z;
+print("H1=",H1);
+H1: hfactor(H1);
+print("C code:");
+cform(H1);
+print("----------------------------------------------------------------");
+
+H2: (y-yc)*z;
+print("H2=",H2);
+H2: hfactor(H2);
+print("C code:");
+cform(H2);
+print("----------------------------------------------------------------");
+
+H3: (x-xc)*(y-yc)*z;
+print("H3=",H3);
+H3: hfactor(H3);
+print("C code:");
+cform(H3);
+print("----------------------------------------------------------------");
+
+H4: z^2;
+print("H4=",H4);
+H4: hfactor(H4);
+print("C code:");
+cform(H4);
+print("----------------------------------------------------------------");
diff --git a/modules/terrain.mod b/modules/terrain.mod
index cc94431..11016c4 100644
--- a/modules/terrain.mod
+++ b/modules/terrain.mod
@@ -120,9 +120,36 @@ static gboolean polygon_contains (Polygon * p, gdouble q[2])
return TRUE;
}
+static gboolean polygon_includes (RSurfaceRect rect, Polygon * p)
+{
+ gdouble q[2];
+ q[0] = rect[0].l; q[1] = rect[1].l;
+ if (!polygon_contains (p, q))
+ return FALSE;
+ q[0] = rect[0].l; q[1] = rect[1].h;
+ if (!polygon_contains (p, q))
+ return FALSE;
+ q[0] = rect[0].h; q[1] = rect[1].l;
+ if (!polygon_contains (p, q))
+ return FALSE;
+ q[0] = rect[0].h; q[1] = rect[1].h;
+ if (!polygon_contains (p, q))
+ return FALSE;
+ return TRUE;
+}
+
+static gboolean polygon_intersects (RSurfaceRect rect, Polygon * p)
+{
+ /* fixme: this could be improved? */
+ return (rect[0].l <= p->max[0] && rect[0].h >= p->min[0] &&
+ rect[1].l <= p->max[1] && rect[1].h >= p->min[1]);
+}
+
typedef struct {
gdouble H[NM+1], m[NM][NM];
gdouble h[NM], he, cond, min, max;
+ GfsRefineTerrain * t;
+ FttCell * cell;
Polygon * p;
gboolean relative;
} RMS;
@@ -136,6 +163,8 @@ static void rms_init (RMS * rms, Polygon * p, gboolean relative)
for (j = 0; j < NM; j++)
rms->m[i][j] = 0.;
rms->p = p;
+ rms->t = p->t;
+ rms->cell = p->cell;
rms->relative = relative;
rms->min = G_MAXDOUBLE;
rms->max = - G_MAXDOUBLE;
@@ -180,6 +209,19 @@ static gdouble cell_value (FttCell * cell, GfsVariable * h[NM], FttVector p)
GFS_VALUE (cell, h[3])*p.x*p.y;
}
+static void cell_coefficients (FttCell * cell, GfsVariable * h[NM], gdouble hp[NM])
+{
+ gdouble size = ftt_cell_size (cell)/2.;
+ FttVector q;
+ ftt_cell_pos (cell, &q);
+ hp[0] = GFS_VALUE (cell, h[0])
+ - (GFS_VALUE (cell, h[1])*q.x + GFS_VALUE (cell, h[2])*q.y
+ - GFS_VALUE (cell, h[3])*q.x*q.y/size)/size;
+ hp[1] = (GFS_VALUE (cell, h[1]) - GFS_VALUE (cell, h[3])*q.y/size)/size;
+ hp[2] = (GFS_VALUE (cell, h[2]) - GFS_VALUE (cell, h[3])*q.x/size)/size;
+ hp[3] = GFS_VALUE (cell, h[3])/(size*size);
+}
+
static void corners_from_parent (FttCell * cell, GfsRefineTerrain * t, gdouble H[4])
{
gdouble size = ftt_cell_size (cell);
@@ -209,7 +251,7 @@ static void variance_check (RMS * rms)
gdouble H[4];
if (rms->relative) {
gdouble H0[4];
- corners_from_parent (rms->p->cell, rms->p->t, H0);
+ corners_from_parent (rms->cell, rms->t, H0);
H[0] = clamp (rms->h[0] + rms->h[1] + rms->h[2] + rms->h[3],
rms->min - H0[0], rms->max - H0[0]);
H[1] = clamp (rms->h[0] - rms->h[1] + rms->h[2] - rms->h[3],
@@ -253,7 +295,7 @@ static void rms_update (RMS * rms)
m[i+NM*j] = rms->m[i][j];
gsl_matrix_view gm = gsl_matrix_view_array (m, NM, NM);
gsl_linalg_SV_decomp_jacobi (&gm.matrix, &gv.matrix, &gs.vector);
- rms->cond = s[NM - 1] > 0. ? s[0]/s[NM - 1] : G_MAXDOUBLE;
+ rms->cond = s[NM - 1] > 0. ? s[0]/s[NM - 1] : G_MAXDOUBLE;
if (rms->cond < 10000.) {
gsl_vector_view gH = gsl_vector_view_array (rms->H, NM);
gsl_vector_view gh = gsl_vector_view_array (rms->h, NM);
@@ -347,19 +389,142 @@ static int rms_add (double p[3], RMS * rms)
return 0;
}
+static void update_terrain_rms (FttCell * cell, GfsRefineTerrain * t, gboolean relative,
+ RMS * rms)
+{
+ Polygon poly;
+ polygon_init (&poly, cell, t);
+ rms_init (rms, &poly, relative);
+ guint i;
+ for (i = 0; i < t->nrs; i++)
+ r_surface_query_region (t->rs[i], poly.min, poly.max, (RSurfaceQuery) rms_add, rms);
+ rms->p = NULL;
+}
+
+static void update_terrain_rms1 (FttCell * cell, GfsRefineTerrain * t, gboolean relative,
+ RMS * rms)
+{
+ Polygon poly;
+ polygon_init (&poly, cell, t);
+ rms_init (rms, &poly, relative);
+ RSurfaceSum s;
+ r_surface_sum_init (&s);
+ guint i;
+ for (i = 0; i < t->nrs; i++)
+ r_surface_query_region_sum (t->rs[i],
+ (RSurfaceCheck) polygon_includes,
+ (RSurfaceCheck) polygon_intersects,
+ &poly, &s);
+ rms->m[0][0] = s.n;
+ if (s.n > 0) {
+ gdouble xc = poly.c.x, yc = poly.c.y;
+ gdouble h = poly.h;
+ /* See terrain.mac for a "maxima" derivation of the terms below */
+ rms->m[0][1] = s.m01 - xc*s.n;
+ rms->m[0][1] /= h;
+ rms->m[0][2] = s.m02 - yc*s.n;
+ rms->m[0][2] /= h;
+ rms->m[0][3] = xc*yc*s.n - s.m01*yc - s.m02*xc + s.m03;
+ rms->m[0][3] /= h*h;
+ rms->m[1][2] = rms->m[0][3];
+ rms->m[1][1] = xc*xc*s.n - 2.*s.m01*xc + s.m11;
+ rms->m[1][1] /= h*h;
+ rms->m[2][2] = yc*yc*s.n - 2.*s.m02*yc + s.m22;
+ rms->m[2][2] /= h*h;
+ rms->m[1][3] = - xc*xc*yc*s.n + 2.*s.m01*xc*yc - s.m11*yc + s.m02*xc*xc - 2.*s.m03*xc + s.m13;
+ rms->m[1][3] /= h*h*h;
+ rms->m[2][3] = - xc*yc*yc*s.n + s.m01*yc*yc + 2.*s.m02*xc*yc - 2.*s.m03*yc - s.m22*xc + s.m23;
+ rms->m[2][3] /= h*h*h;
+ rms->m[3][3] = xc*xc*yc*yc*s.n - 2.*s.m01*xc*yc*yc + s.m11*yc*yc - 2.*s.m02*xc*xc*yc
+ + 4.*s.m03*xc*yc - 2.*s.m13*yc + s.m22*xc*xc - 2.*s.m23*xc + s.m33;
+ rms->m[3][3] /= h*h*h*h;
+ if (rms->relative) {
+ double hp[NM];
+ cell_coefficients (ftt_cell_parent (rms->cell), rms->t->h, hp);
+ rms->H[0] = s.H0 - s.n*hp[0] - s.m01*hp[1] - s.m02*hp[2] - s.m03*hp[3];
+ /* See terrain.mac for a "maxima" derivation of the terms below */
+ rms->H[1] = (s.H1 - xc*s.H0
+ + s.m03*(hp[3]*xc - hp[2])
+ + s.m01*(hp[1]*xc - hp[0])
+ + hp[2]*s.m02*xc
+ + hp[0]*xc*s.n
+ - hp[3]*s.m13
+ - hp[1]*s.m11);
+ rms->H[2] = (s.H2 - yc*s.H0
+ + s.m03*(hp[3]*yc - hp[1])
+ + s.m02*(hp[2]*yc - hp[0])
+ + hp[1]*s.m01*yc
+ + hp[0]*yc*s.n
+ - hp[3]*s.m23
+ - hp[2]*s.m22);
+ rms->H[3] = (s.H3 - xc*s.H2 - yc*s.H1 + xc*yc*s.H0
+ - s.m03*(hp[3]*xc*yc - hp[2]*yc - hp[1]*xc + hp[0])
+ + s.m13*(hp[3]*yc - hp[1])
+ - s.m02*xc*(hp[2]*yc - hp[0])
+ - s.m01*(hp[1]* xc - hp[0])*yc
+ - hp[0]*xc*yc*s.n
+ + hp[1]*s.m11*yc
+ + s.m23*(hp[3]*xc - hp[2])
+ + hp[2]*s.m22*xc
+ - hp[3]*s.m33);
+ rms->H[4] = (s.H4 - 2.*hp[3]*s.H3 - 2.*hp[2]*s.H2 - 2.*hp[1]*s.H1 - 2.*hp[0]*s.H0
+ + hp[3]*hp[3]*s.m33
+ + 2.*hp[2]*hp[3]*s.m23
+ + hp[2]*hp[2]*s.m22
+ + 2.*hp[1]*hp[3]*s.m13
+ + hp[1]*hp[1]*s.m11
+ + 2.*(hp[0]*hp[3] + hp[1]*hp[2])*s.m03
+ + 2.*hp[0]*hp[2]*s.m02
+ + 2.*hp[0]*hp[1]*s.m01
+ + hp[0]*hp[0]*s.n);
+ }
+ else {
+ rms->H[0] = s.H0;
+ rms->H[1] = s.H1 - xc*s.H0;
+ rms->H[2] = s.H2 - yc*s.H0;
+ rms->H[3] = s.H3- xc*s.H2 - yc*s.H1 + xc*yc*s.H0;
+ rms->H[4] = s.H4;
+ }
+ rms->H[1] /= h;
+ rms->H[2] /= h;
+ rms->H[3] /= h*h;
+ rms->max = s.Hmax;
+ rms->min = s.Hmin;
+ }
+ rms->p = NULL;
+}
+
static void update_terrain (FttCell * cell, GfsRefineTerrain * t)
{
- g_assert (GFS_VALUE (cell, t->type) == REFINED);
- Polygon p;
- polygon_init (&p, cell, t);
RMS rms;
- rms_init (&rms, &p, ftt_cell_parent (cell) != NULL);
guint i;
- for (i = 0; i < t->nrs; i++)
- r_surface_query_region (t->rs[i], p.min, p.max, (RSurfaceQuery) rms_add, &rms);
+ g_assert (GFS_VALUE (cell, t->type) == REFINED);
+ update_terrain_rms (cell, t, ftt_cell_parent (cell) != NULL, &rms);
+#if 1
+ RMS rms1;
+ update_terrain_rms1 (cell, t, ftt_cell_parent (cell) != NULL, &rms1);
+ {
+ guint i, j;
+ for (i = 0; i < NM; i++)
+ for (j = 0; j <= i; j++)
+ fprintf (stderr, "m%d%d: %g %g\n", j, i, rms.m[j][i], rms1.m[j][i]);
+ }
+ {
+ guint i;
+ for (i = 0; i < NM + 1; i++)
+ fprintf (stderr, "H%s%d: %g %g\n", rms.relative ? "*" : "", i, rms.H[i], rms1.H[i]);
+ }
+ if (rms1.m[0][0] > 0.) {
+ fprintf (stderr, "min: %g %g\n", rms.min, rms1.min);
+ fprintf (stderr, "max: %g %g\n", rms.max, rms1.max);
+ }
+ rms_update (&rms1);
+#endif
rms_update (&rms);
- for (i = 0; i < NM; i++)
+ for (i = 0; i < NM; i++) {
GFS_VALUE (cell, t->h[i]) = rms.h[i];
+ fprintf (stderr, "h%d: %g %g\n", i, rms.h[i], rms1.h[i]);
+ }
GFS_VALUE (cell, t->he) = rms.he;
GFS_VALUE (cell, t->hn) = rms.m[0][0];
GFS_VALUE (cell, t->type) = RAW;
@@ -420,15 +585,18 @@ static gdouble corner_value (GfsRefineTerrain * t, FttVector * p, gdouble eps, g
static void update_error_estimate (FttCell * cell, GfsRefineTerrain * t, gboolean relative)
{
if (GFS_VALUE (cell, t->hn) > 0.) {
- Polygon poly;
- polygon_init (&poly, cell, t);
RMS rms;
- rms_init (&rms, &poly, relative);
guint i;
- for (i = 0; i < t->nrs; i++)
- r_surface_query_region (t->rs[i], poly.min, poly.max, (RSurfaceQuery) rms_add, &rms);
+ update_terrain_rms (cell, t, relative, &rms);
for (i = 0; i < NM; i++)
rms.h[i] = GFS_VALUE (cell, t->h[i]);
+#if 1
+ RMS rms1;
+ update_terrain_rms1 (cell, t, relative, &rms1);
+ for (i = 0; i < NM; i++)
+ rms1.h[i] = GFS_VALUE (cell, t->h[i]);
+ fprintf (stderr, "he: %g %g\n", rms_minimum (&rms), rms_minimum (&rms1));
+#endif
GFS_VALUE (cell, t->he) = rms_minimum (&rms);
}
else
--
Gerris Flow Solver
More information about the debian-science-commits
mailing list