[SCM] Gerris Flow Solver branch, upstream, updated. b3aa46814a06c9cb2912790b23916ffb44f1f203
Stephane Popinet
popinet at users.sf.net
Fri May 15 02:54:10 UTC 2009
The following commit has been merged in the upstream branch:
commit ccd172d5a48779ae7bdb6424aad97e9b05dcf8fc
Author: Stephane Popinet <popinet at users.sf.net>
Date: Tue Mar 27 10:23:04 2007 +1000
Adaptivity along solid boundaries should now work (but not for 2D3 yet)
darcs-hash:20070327002304-d4795-cde68d37175d2c2ac1b516717bb62eedab00cb71.gz
diff --git a/src/adaptive.c b/src/adaptive.c
index ca485cd..7b9e627 100644
--- a/src/adaptive.c
+++ b/src/adaptive.c
@@ -69,8 +69,8 @@ void gfs_cell_fine_init (FttCell * parent, GfsDomain * domain)
gfs_cell_init (parent, domain);
- /* fixme: refinement of mixed cell is not implemented (yet) */
- g_assert (GFS_CELL_IS_BOUNDARY (parent) || GFS_IS_FLUID (parent));
+ if (!GFS_CELL_IS_BOUNDARY (parent) && GFS_IS_MIXED (parent))
+ gfs_solid_coarse_fine (parent);
i = domain->variables;
while (i) {
diff --git a/src/fluid.c b/src/fluid.c
index 168284e..514be9f 100644
--- a/src/fluid.c
+++ b/src/fluid.c
@@ -1554,7 +1554,8 @@ void gfs_cell_coarse_fine (FttCell * parent, GfsVariable * v)
ftt_cell_children (parent, &child);
for (n = 0; n < FTT_CELLS; n++)
- GFS_VARIABLE (child.c[n], v->i) = GFS_VARIABLE (parent, v->i);
+ if (child.c[n])
+ GFS_VARIABLE (child.c[n], v->i) = GFS_VARIABLE (parent, v->i);
if (!GFS_CELL_IS_BOUNDARY (parent)) {
FttVector g;
@@ -1563,13 +1564,14 @@ void gfs_cell_coarse_fine (FttCell * parent, GfsVariable * v)
for (c = 0; c < FTT_DIMENSION; c++)
(&g.x)[c] = gfs_center_van_leer_gradient (parent, c, v->i);
- for (n = 0; n < FTT_CELLS; n++) {
- FttVector p;
-
- ftt_cell_relative_pos (child.c[n], &p);
- for (c = 0; c < FTT_DIMENSION; c++)
- GFS_VARIABLE (child.c[n], v->i) += (&p.x)[c]*(&g.x)[c];
- }
+ for (n = 0; n < FTT_CELLS; n++)
+ if (child.c[n]) {
+ FttVector p;
+
+ ftt_cell_relative_pos (child.c[n], &p);
+ for (c = 0; c < FTT_DIMENSION; c++)
+ GFS_VARIABLE (child.c[n], v->i) += (&p.x)[c]*(&g.x)[c];
+ }
}
}
diff --git a/src/solid.c b/src/solid.c
index 38b6507..50220e7 100644
--- a/src/solid.c
+++ b/src/solid.c
@@ -1305,3 +1305,161 @@ void gfs_face_ca (const FttCellFace * face, FttVector * ca)
#endif /* 3D */
}
}
+
+#if !FTT_2D
+static void outer_fractions_coarse_fine (FttCell * parent, FttDirection d)
+{
+ GfsSolidVector * solid = GFS_STATE (parent)->solid;
+ FttComponent c1 = d < 4 ? 2 : 0, c2 = d < 2 || d > 3 ? 1 : 0;
+ gdouble nm;
+ FttVector m;
+
+ m.x = solid->s[2*c1 + 1] - solid->s[2*c1]; nm = fabs (m.x);
+ m.y = solid->s[2*c2 + 1] - solid->s[2*c2]; nm += fabs (m.y);
+ if (nm > 0.) {
+ m.x /= nm;
+ m.y /= nm;
+ }
+ else
+ m.x = 1.;
+ gdouble alpha = gfs_line_alpha (&m, solid->s[d]);
+ gdouble ss = 0.;
+
+ FttCellChildren child;
+ guint i, n = ftt_cell_children_direction (parent, d, &child);
+ for (i = 0; i < n; i++)
+ if (child.c[i]) {
+ if (GFS_IS_MIXED (child.c[i])) {
+ GfsSolidVector * s = GFS_STATE (child.c[i])->solid;
+ gdouble alpha1 = alpha;
+ FttVector p;
+
+ ftt_cell_relative_pos (child.c[i], &p);
+ alpha1 -= m.x*(0.25 + (&p.x)[c1]);
+ alpha1 -= m.y*(0.25 + (&p.x)[c2]);
+
+ s->s[d] = gfs_line_area (&m, 2.*alpha1);
+ ss += s->s[d];
+ }
+ else
+ ss += 1.;
+ }
+ /* fixme: this should not happen
+ * It happens in configurations where children cells are not cut by
+ * the VOF approximation but should have non-zero surface
+ * fractions */
+ if (fabs (solid->s[d] - ss/n) > 1e-5)
+ g_warning ("inconsistent surface fractions %d %f %f %f\n", d, solid->s[d], ss/n,
+ fabs (solid->s[d] - ss/n));
+}
+#endif /* 3D */
+
+/**
+ * gfs_solid_coarse_fine:
+ * @parent: a mixed #FttCell with children.
+ *
+ * Fills the solid properties of the children of @parent.
+ * Destroys all children entirely contained in the solid.
+ */
+void gfs_solid_coarse_fine (FttCell * parent)
+{
+#if FTT_2D3
+ g_assert_not_implemented ();
+#endif
+ g_return_if_fail (parent);
+ g_return_if_fail (GFS_IS_MIXED (parent));
+ g_return_if_fail (!FTT_CELL_IS_LEAF (parent));
+
+ GfsSolidVector * solid = GFS_STATE (parent)->solid;
+ FttVector m;
+ FttComponent c;
+ gdouble n = 0;
+
+ for (c = 0; c < FTT_DIMENSION; c++) {
+ (&m.x)[c] = solid->s[2*c + 1] - solid->s[2*c];
+ n += fabs ((&m.x)[c]);
+ }
+ if (n > 0.)
+ for (c = 0; c < FTT_DIMENSION; c++)
+ (&m.x)[c] /= n;
+ else
+ m.x = 1.;
+ gdouble alpha = gfs_plane_alpha (&m, solid->a);
+
+ gdouble h = ftt_cell_size (parent)/2.;
+ guint level = ftt_cell_level (parent) + 1;
+ FttCellChildren child;
+ guint i;
+ ftt_cell_children (parent, &child);
+ for (i = 0; i < FTT_CELLS; i++) {
+ gdouble alpha1 = alpha;
+ FttVector p;
+
+ ftt_cell_relative_pos (child.c[i], &p);
+ for (c = 0; c < FTT_DIMENSION; c++)
+ alpha1 -= (&m.x)[c]*(0.25 + (&p.x)[c]);
+
+ if (GFS_STATE (child.c[i])->solid) {
+ g_free (GFS_STATE (child.c[i])->solid);
+ GFS_STATE (child.c[i])->solid = NULL;
+ }
+
+ gdouble a = gfs_plane_volume (&m, 2.*alpha1);
+ if (a > 0. && a < 1.) {
+ GfsSolidVector * s = GFS_STATE (child.c[i])->solid = g_malloc (sizeof (GfsSolidVector));
+ s->a = a;
+
+ ftt_cell_pos (child.c[i], &p);
+ gfs_plane_center (&m, 2.*alpha1, a, &s->cm);
+ for (c = 0; c < FTT_DIMENSION; c++)
+ (&s->cm.x)[c] = (&p.x)[c] + h*((&s->cm.x)[c] - 0.5);
+ g_assert (gfs_vof_plane_center (child.c[i], &m, 2.*alpha1, &s->ca));
+
+ FttDirection d;
+ FttCellNeighbors n;
+ ftt_cell_neighbors (child.c[i], &n);
+ for (d = 0; d < FTT_NEIGHBORS; d++)
+ if (!n.c[d])
+ s->s[d] = 0.;
+ else if (GFS_IS_MIXED (n.c[d]) && ftt_cell_level (n.c[d]) == level)
+ s->s[d] = GFS_STATE (n.c[d])->solid->s[FTT_OPPOSITE_DIRECTION (d)];
+ else if (!ftt_cell_neighbor_is_brother (child.c[i], d) && GFS_IS_FLUID (n.c[d]))
+ s->s[d] = 1.;
+ else {
+#if FTT_2D
+ gdouble f;
+ FttComponent c1 = d > 1, c2 = !c1;
+
+ if ((&m.x)[c2] == 0.) f = 0.;
+ else {
+ f = (2.*alpha1 - (&m.x)[c1]*!(d % 2))/(&m.x)[c2];
+ if (f < 0.) f = 0.; else if (f > 1.) f = 1.;
+ if ((&m.x)[c2] < 0.)
+ f = 1. - f;
+ }
+ s->s[d] = f;
+#else /* 3D */
+ /* only initialises "inner" fractions */
+ if (ftt_cell_neighbor_is_brother (child.c[i], d)) {
+ FttComponent c1 = (d/2 + 1) % 3, c2 = (d/2 + 2) % 3;
+ FttVector mp;
+ mp.x = (&m.x)[c1];
+ mp.y = (&m.x)[c2];
+ s->s[d] = gfs_line_area (&mp, d % 2 ? 2.*alpha1 : 2.*alpha1 - (&m.x)[d/2]);
+ }
+#endif /* 3D */
+ }
+ }
+ else if (a == 0.)
+ ftt_cell_destroy (child.c[i], (FttCellCleanupFunc) gfs_cell_cleanup, NULL);
+ }
+
+#if !FTT_2D
+ FttCellNeighbors neighbor;
+ FttDirection d;
+ ftt_cell_neighbors (parent, &neighbor);
+ for (d = 0; d < FTT_NEIGHBORS; d++)
+ if (neighbor.c[d] && FTT_CELL_IS_LEAF (neighbor.c[d]) && !GFS_IS_FLUID (neighbor.c[d]))
+ outer_fractions_coarse_fine (parent, d);
+#endif /* 3D */
+}
diff --git a/src/solid.h b/src/solid.h
index 1cfda72..5230d06 100644
--- a/src/solid.h
+++ b/src/solid.h
@@ -50,6 +50,7 @@ void gfs_solid_normal (const FttCell * cell,
FttVector * n);
void gfs_face_ca (const FttCellFace * face,
FttVector * ca);
+void gfs_solid_coarse_fine (FttCell * parent);
#ifdef __cplusplus
}
diff --git a/src/tension.c b/src/tension.c
index 0a67fd4..6fc778e 100644
--- a/src/tension.c
+++ b/src/tension.c
@@ -566,14 +566,15 @@ static void curvature_coarse_fine (FttCell * parent, GfsVariable * v)
guint n;
ftt_cell_children (parent, &child);
- for (n = 0; n < FTT_CELLS; n++) {
- GfsVariableCurvature * k = GFS_VARIABLE_CURVATURE (v);
- gdouble f = GFS_VARIABLE (child.c[n], k->f->i);
-
- GFS_VARIABLE (child.c[n], v->i) = GFS_VARIABLE (parent, v->i);
- if (!GFS_IS_FULL (f))
- g_assert (fabs (GFS_VARIABLE (child.c[n], v->i)) < G_MAXDOUBLE);
- }
+ for (n = 0; n < FTT_CELLS; n++)
+ if (child.c[n]) {
+ GfsVariableCurvature * k = GFS_VARIABLE_CURVATURE (v);
+ gdouble f = GFS_VARIABLE (child.c[n], k->f->i);
+
+ GFS_VARIABLE (child.c[n], v->i) = GFS_VARIABLE (parent, v->i);
+ if (!GFS_IS_FULL (f))
+ g_assert (fabs (GFS_VARIABLE (child.c[n], v->i)) < G_MAXDOUBLE);
+ }
}
static void curvature_fine_coarse (FttCell * parent, GfsVariable * v)
diff --git a/src/vof.c b/src/vof.c
index 0d8cd27..1bcae33 100644
--- a/src/vof.c
+++ b/src/vof.c
@@ -740,6 +740,8 @@ static void vof_coarse_fine (FttCell * parent, GfsVariable * v)
if (GFS_IS_FULL (f))
for (i = 0; i < FTT_CELLS; i++) {
FttComponent c;
+ if (!child.c[i])
+ g_assert_not_implemented ();
GFS_VARIABLE (child.c[i], v->i) = f;
for (c = 1; c < FTT_DIMENSION; c++)
GFS_VARIABLE (child.c[i], t->m[c]->i) = 0.;
@@ -757,6 +759,8 @@ static void vof_coarse_fine (FttCell * parent, GfsVariable * v)
FttComponent c;
FttVector p;
+ if (!child.c[i])
+ g_assert_not_implemented ();
ftt_cell_relative_pos (child.c[i], &p);
for (c = 0; c < FTT_DIMENSION; c++) {
alpha1 -= (&m.x)[c]*(0.25 + (&p.x)[c]);
--
Gerris Flow Solver
More information about the debian-science-commits
mailing list