[SCM] Gerris Flow Solver branch, upstream, updated. e8f73a07832050124d2b8bf6c6f35b33180e65a8
Stephane Popinet
popinet at users.sf.net
Tue Nov 24 12:25:10 UTC 2009
The following commit has been merged in the upstream branch:
commit 79a35c8947e9919b97ed3efa378bec03179391ed
Author: Stephane Popinet <popinet at users.sf.net>
Date: Wed Oct 21 10:13:04 2009 +1100
Support for "rotated" periodic boundaries
darcs-hash:20091020231304-d4795-943b51d73f824928bd610480814c5eb8efa39f22.gz
diff --git a/src/boundary.c b/src/boundary.c
index b7cb5fd..26306b0 100644
--- a/src/boundary.c
+++ b/src/boundary.c
@@ -1164,7 +1164,8 @@ static void send (GfsBoundary * bb)
{
GfsBoundaryPeriodic * boundary = GFS_BOUNDARY_PERIODIC (bb);
g_assert (boundary->matching);
- GfsBoundaryPeriodic * matching = GFS_BOUNDARY_PERIODIC (boundary->matching->neighbor[bb->d]);
+ GfsBoundaryPeriodic * matching =
+ GFS_BOUNDARY_PERIODIC (boundary->matching->neighbor[boundary->d]);
g_assert (GFS_IS_BOUNDARY_PERIODIC (matching));
g_assert (boundary->sndcount <= boundary->sndbuf->len);
@@ -1328,6 +1329,8 @@ static void gfs_boundary_periodic_init (GfsBoundaryPeriodic * boundary)
boundary->sndbuf = g_array_new (FALSE, FALSE, sizeof (gdouble));
boundary->rcvbuf = g_array_new (FALSE, FALSE, sizeof (gdouble));
boundary->sndcount = boundary->rcvcount = 0;
+
+ boundary->rotate = 0.;
}
GfsBoundaryClass * gfs_boundary_periodic_class (void)
@@ -1351,6 +1354,15 @@ GfsBoundaryClass * gfs_boundary_periodic_class (void)
return klass;
}
+/**
+ * gfs_boundary_periodic_new:
+ * @klass: a #GfsBoundaryClass.
+ * @box: a #GfsBox.
+ * @d: a #FttDirection.
+ * @matching: a #GfsBox.
+ *
+ * Returns: a new #GfsBoundaryPeriodic connecting @box with @matching in direction @d.
+ */
GfsBoundaryPeriodic * gfs_boundary_periodic_new (GfsBoundaryClass * klass,
GfsBox * box,
FttDirection d,
@@ -1361,6 +1373,58 @@ GfsBoundaryPeriodic * gfs_boundary_periodic_new (GfsBoundaryClass * klass,
boundary = GFS_BOUNDARY_PERIODIC (gfs_boundary_new (klass, box, d));
set_buffers_size (boundary);
boundary->matching = matching;
+ boundary->d = FTT_OPPOSITE_DIRECTION (d);
+
+ return boundary;
+}
+
+static void center_periodic_rotate (FttCellFace * face, GfsBc * b)
+{
+ GfsBoundaryPeriodic * boundary_periodic = GFS_BOUNDARY_PERIODIC (b->b);
+
+ g_assert (boundary_periodic->sndcount < boundary_periodic->sndbuf->len);
+ g_assert (ftt_face_type (face) == FTT_FINE_FINE);
+ g_assert (!FTT_CELL_IS_LEAF (face->cell) || FTT_CELL_IS_LEAF (face->neighbor));
+
+ if (b->v->component < 2) { /* 2D-vector-rotation only */
+ FttComponent c = FTT_ORTHOGONAL_COMPONENT (b->v->component);
+ g_assert (b->v->vector[c]);
+ g_array_index (boundary_periodic->sndbuf, gdouble, boundary_periodic->sndcount++) =
+ (2.*c - 1.)*boundary_periodic->rotate*GFS_VALUE (face->neighbor, b->v->vector[c]);
+ }
+ else
+ g_array_index (boundary_periodic->sndbuf, gdouble, boundary_periodic->sndcount++) =
+ GFS_VALUE (face->neighbor, b->v);
+}
+
+/**
+ * gfs_boundary_periodic_rotate_new:
+ * @klass: a #GfsBoundaryClass.
+ * @box: a #GfsBox.
+ * @d: a #FttDirection.
+ * @matching: a #GfsBox.
+ * @rotate: a #FttDirection.
+ * @orientation: the orientation (+1 or -1).
+ *
+ * Returns: a new "rotated" #GfsBoundaryPeriodic connecting @box in
+ * direction @d with @matching in direction @rotate, oriented using
+ * @orientation.
+ */
+GfsBoundaryPeriodic * gfs_boundary_periodic_rotate_new (GfsBoundaryClass * klass,
+ GfsBox * box,
+ FttDirection d,
+ GfsBox * matching,
+ FttDirection rotate,
+ gdouble orientation)
+{
+ GfsBoundaryPeriodic * boundary;
+
+ boundary = gfs_boundary_periodic_new (klass, box, d, matching);
+ boundary->d = rotate;
+ boundary->rotate = orientation;
+
+ GfsBc * b = GFS_BOUNDARY (boundary)->default_bc;
+ b->bc = b->homogeneous_bc = (FttFaceTraverseFunc) center_periodic_rotate;
return boundary;
}
@@ -1370,6 +1434,8 @@ GfsBoundaryPeriodic * gfs_boundary_periodic_new (GfsBoundaryClass * klass,
static void gfs_gedge_write (GtsObject * object, FILE * fp)
{
fprintf (fp, " %s", ftt_direction_name [GFS_GEDGE (object)->d]);
+ if (GFS_GEDGE (object)->rotate < FTT_NEIGHBORS)
+ fprintf (fp, " %s", ftt_direction_name [GFS_GEDGE (object)->rotate]);
}
static void gfs_gedge_read (GtsObject ** o, GtsFile * fp)
@@ -1387,6 +1453,15 @@ static void gfs_gedge_read (GtsObject ** o, GtsFile * fp)
return;
}
gts_file_next_token (fp);
+ if (fp->type == GTS_STRING) {
+ e->rotate = ftt_direction_from_name (fp->token->str);
+ if (e->rotate >= FTT_NEIGHBORS) {
+ gts_file_error (fp, "unknown direction `%s'", fp->token->str);
+ e->rotate = -1;
+ return;
+ }
+ gts_file_next_token (fp);
+ }
}
static void gfs_gedge_class_init (GtsObjectClass * klass)
@@ -1397,7 +1472,7 @@ static void gfs_gedge_class_init (GtsObjectClass * klass)
static void gfs_gedge_init (GfsGEdge * object)
{
- object->d = -1;
+ object->d = object->rotate = FTT_NEIGHBORS;
}
GfsGEdgeClass * gfs_gedge_class (void)
@@ -1442,22 +1517,32 @@ void gfs_gedge_link_boxes (GfsGEdge * edge)
b2 = GFS_BOX (GTS_GEDGE (edge)->n2);
g_return_if_fail (b1->neighbor[edge->d] == NULL);
- g_return_if_fail (b2->neighbor[FTT_OPPOSITE_DIRECTION (edge->d)] == NULL);
-
- GtsObject * periodic = GTS_OBJECT (b1);
- while (periodic && GFS_IS_BOX (periodic) && GFS_BOX (periodic) != b2)
- periodic = GFS_BOX (periodic)->neighbor[FTT_OPPOSITE_DIRECTION (edge->d)];
-
- if (GFS_BOX (periodic) == b2) {
- gfs_boundary_periodic_new (gfs_boundary_periodic_class (), b1, edge->d, b2);
- gfs_boundary_periodic_new (gfs_boundary_periodic_class (), b2,
- FTT_OPPOSITE_DIRECTION (edge->d), b1);
+
+ if (edge->rotate < FTT_NEIGHBORS) {
+ g_return_if_fail (b2->neighbor[edge->rotate] == NULL);
+ gfs_boundary_periodic_rotate_new (gfs_boundary_periodic_class (),
+ b1, edge->d, b2, edge->rotate, 1.);
+ gfs_boundary_periodic_rotate_new (gfs_boundary_periodic_class (),
+ b2, edge->rotate, b1, edge->d, -1.);
}
else {
- ftt_cell_set_neighbor (b1->root, b2->root, edge->d,
- (FttCellInitFunc) gfs_cell_init, gfs_box_domain (b1));
- b1->neighbor[edge->d] = GTS_OBJECT (b2);
- b2->neighbor[FTT_OPPOSITE_DIRECTION (edge->d)] = GTS_OBJECT (b1);
+ g_return_if_fail (b2->neighbor[FTT_OPPOSITE_DIRECTION (edge->d)] == NULL);
+
+ GtsObject * periodic = GTS_OBJECT (b1);
+ while (periodic && GFS_IS_BOX (periodic) && GFS_BOX (periodic) != b2)
+ periodic = GFS_BOX (periodic)->neighbor[FTT_OPPOSITE_DIRECTION (edge->d)];
+
+ if (GFS_BOX (periodic) == b2) {
+ gfs_boundary_periodic_new (gfs_boundary_periodic_class (), b1, edge->d, b2);
+ gfs_boundary_periodic_new (gfs_boundary_periodic_class (), b2,
+ FTT_OPPOSITE_DIRECTION (edge->d), b1);
+ }
+ else {
+ ftt_cell_set_neighbor (b1->root, b2->root, edge->d,
+ (FttCellInitFunc) gfs_cell_init, gfs_box_domain (b1));
+ b1->neighbor[edge->d] = GTS_OBJECT (b2);
+ b2->neighbor[FTT_OPPOSITE_DIRECTION (edge->d)] = GTS_OBJECT (b1);
+ }
}
}
diff --git a/src/boundary.h b/src/boundary.h
index 677bed1..217ca7f 100644
--- a/src/boundary.h
+++ b/src/boundary.h
@@ -247,8 +247,11 @@ struct _GfsBoundaryPeriodic {
GfsBoundary parent;
GfsBox * matching;
+ FttDirection d;
GArray * sndbuf, * rcvbuf;
guint sndcount, rcvcount;
+
+ gdouble rotate;
};
#define GFS_BOUNDARY_PERIODIC(obj) GTS_OBJECT_CAST (obj,\
@@ -262,6 +265,12 @@ GfsBoundaryPeriodic * gfs_boundary_periodic_new (GfsBoundaryClass * klass,
GfsBox * box,
FttDirection d,
GfsBox * matching);
+GfsBoundaryPeriodic * gfs_boundary_periodic_rotate_new (GfsBoundaryClass * klass,
+ GfsBox * box,
+ FttDirection d,
+ GfsBox * matching,
+ FttDirection rotate,
+ gdouble orientation);
/* GfsGEdge: Header */
@@ -271,7 +280,7 @@ typedef struct _GfsGEdgeClass GfsGEdgeClass;
struct _GfsGEdge {
GtsGEdge parent;
- FttDirection d;
+ FttDirection d, rotate;
};
struct _GfsGEdgeClass {
diff --git a/src/event.c b/src/event.c
index 1b18030..7d45a79 100644
--- a/src/event.c
+++ b/src/event.c
@@ -624,9 +624,16 @@ static gboolean gfs_init_event (GfsEvent * event, GfsSimulation * sim)
while (i) {
VarFunc * vf = i->data;
- gfs_traverse_and_bc (GFS_DOMAIN (sim), FTT_PRE_ORDER, FTT_TRAVERSE_LEAFS, -1,
- (FttCellTraverseFunc) init_vf, vf,
- vf->v, vf->v);
+ gfs_domain_cell_traverse (GFS_DOMAIN (sim), FTT_PRE_ORDER, FTT_TRAVERSE_LEAFS, -1,
+ (FttCellTraverseFunc) init_vf, vf);
+ i = i->next;
+ }
+ /* boundary conditions need to be called in a separate loop so
+ that they are correctly applied for vector quantities */
+ i = GFS_INIT (event)->f;
+ while (i) {
+ VarFunc * vf = i->data;
+ gfs_domain_bc (GFS_DOMAIN (sim), FTT_TRAVERSE_LEAFS, -1, vf->v);
i = i->next;
}
return TRUE;
diff --git a/src/graphic.c b/src/graphic.c
index 2590185..ed3c354 100644
--- a/src/graphic.c
+++ b/src/graphic.c
@@ -1718,10 +1718,9 @@ static GList * grow_curve (GfsDomain * domain,
gpointer data[2];
path_class = GTS_POINT_CLASS (gfs_twisted_vertex_class ());
- for (c = 0; c < FTT_DIMENSION; c++) {
+ for (c = 0; c < FTT_DIMENSION; c++)
vort[c] = gfs_temporary_variable (domain);
- gfs_variable_set_vector (vort[c], c);
- }
+ gfs_variable_set_vector (vort, FTT_DIMENSION);
data[0] = vort;
data[1] = U;
gfs_domain_cell_traverse (domain, FTT_PRE_ORDER, FTT_TRAVERSE_LEAFS, -1,
diff --git a/src/moving.c b/src/moving.c
index a0c2ba4..4591891 100644
--- a/src/moving.c
+++ b/src/moving.c
@@ -822,14 +822,13 @@ static void simulation_moving_run (GfsSimulation * sim)
FttComponent c;
for (c = 0; c < FTT_DIMENSION; c++) {
gmac[c] = gfs_temporary_variable (domain);
- gfs_variable_set_vector (gmac[c], c);
- if (sim->advection_params.gc) {
+ if (sim->advection_params.gc)
g[c] = gfs_temporary_variable (domain);
- gfs_variable_set_vector (g[c], c);
- }
else
g[c] = gmac[c];
}
+ gfs_variable_set_vector (gmac, FTT_DIMENSION);
+ gfs_variable_set_vector (g, FTT_DIMENSION);
gfs_simulation_refine (sim);
gfs_simulation_init (sim);
diff --git a/src/ocean.c b/src/ocean.c
index eb2081e..8f27014 100644
--- a/src/ocean.c
+++ b/src/ocean.c
@@ -113,10 +113,9 @@ static void gfs_correct_normal_velocities_weighted (GfsDomain * domain,
g_return_if_fail (p != NULL);
g_return_if_fail (g != NULL);
- for (c = 0; c < dimension; c++) {
+ for (c = 0; c < dimension; c++)
g[c] = gfs_temporary_variable (domain);
- gfs_variable_set_vector (g[c], c);
- }
+ gfs_variable_set_vector (g, dimension);
data[0] = g;
data[1] = &dimension;
gfs_domain_cell_traverse (domain, FTT_PRE_ORDER, FTT_TRAVERSE_LEAFS, -1,
diff --git a/src/river.c b/src/river.c
index f463320..e2bb06e 100644
--- a/src/river.c
+++ b/src/river.c
@@ -462,9 +462,8 @@ static void river_init (GfsRiver * r)
r->v1[0] = gfs_domain_add_variable (domain, NULL, NULL);
r->v1[1] = gfs_domain_add_variable (domain, NULL, NULL);
- r->v1[1]->component = FTT_X;
r->v1[2] = gfs_domain_add_variable (domain, NULL, NULL);
- r->v1[2]->component = FTT_Y;
+ gfs_variable_set_vector (&r->v1[1], 2);
r->dv[0][0] = gfs_domain_add_variable (domain, "Px", "x-component of the thickness gradient");
r->dv[1][0] = gfs_domain_add_variable (domain, "Py", "y-component of the thickness gradien");
diff --git a/src/simulation.c b/src/simulation.c
index ca557e2..fc43f58 100644
--- a/src/simulation.c
+++ b/src/simulation.c
@@ -407,14 +407,13 @@ static void simulation_run (GfsSimulation * sim)
FttComponent c;
for (c = 0; c < FTT_DIMENSION; c++) {
gmac[c] = gfs_temporary_variable (domain);
- gfs_variable_set_vector (gmac[c], c);
- if (sim->advection_params.gc) {
+ if (sim->advection_params.gc)
g[c] = gfs_temporary_variable (domain);
- gfs_variable_set_vector (g[c], c);
- }
else
g[c] = gmac[c];
}
+ gfs_variable_set_vector (gmac, FTT_DIMENSION);
+ gfs_variable_set_vector (g, FTT_DIMENSION);
gfs_simulation_refine (sim);
gfs_simulation_init (sim);
@@ -862,17 +861,16 @@ static void simulation_init (GfsSimulation * object)
v = gfs_domain_add_variable (domain, "Pmac", "MAC projection pressure");
v->centered = TRUE;
v->units = 2.;
- v = gfs_domain_add_variable (domain, "U", "x-component of the velocity");
- gfs_variable_set_vector (v, FTT_X);
- v->units = 1.;
- v = gfs_domain_add_variable (domain, "V", "y-component of the velocity");
- gfs_variable_set_vector (v, FTT_Y);
- v->units = 1.;
+ GfsVariable * u[FTT_DIMENSION];
+ u[0] = gfs_domain_add_variable (domain, "U", "x-component of the velocity");
+ u[0]->units = 1.;
+ u[1] = gfs_domain_add_variable (domain, "V", "y-component of the velocity");
+ u[1]->units = 1.;
#if (!FTT_2D)
- v = gfs_domain_add_variable (domain, "W", "z-component of the velocity");
- gfs_variable_set_vector (v, FTT_Z);
- v->units = 1.;
+ u[2] = gfs_domain_add_variable (domain, "W", "z-component of the velocity");
+ u[2]->units = 1.;
#endif /* FTT_3D */
+ gfs_variable_set_vector (u, FTT_DIMENSION);
GfsDerivedVariableInfo * dv = derived_variable;
while (dv->name) {
diff --git a/src/tension.c b/src/tension.c
index 75bd505..44e38ac 100644
--- a/src/tension.c
+++ b/src/tension.c
@@ -670,10 +670,9 @@ static void variable_curvature_from_distance (GfsEvent * event, GfsSimulation *
if (GFS_IS_AXI (sim))
g_assert_not_implemented ();
- for (c = 0; c < FTT_DIMENSION + 1; c++) {
+ for (c = 0; c < FTT_DIMENSION + 1; c++)
n[c] = gfs_temporary_variable (domain);
- gfs_variable_set_vector (n[c], c);
- }
+ gfs_variable_set_vector (n, FTT_DIMENSION);
data[0] = n;
data[1] = event;
gfs_domain_cell_traverse (domain, FTT_PRE_ORDER, FTT_TRAVERSE_LEAFS, -1,
diff --git a/src/variable.c b/src/variable.c
index 1d1f240..00cef75 100644
--- a/src/variable.c
+++ b/src/variable.c
@@ -233,6 +233,28 @@ void gfs_variables_swap (GfsVariable * v1, GfsVariable * v2)
i = v1->i; v1->i = v2->i; v2->i = i;
}
+/**
+ * gfs_variable_set_vector:
+ * @v: the components of the vector.
+ * @n: the vector dimension.
+ *
+ * Sets @v[0],..., at v[n-1] as components of a vector quantity.
+ */
+void gfs_variable_set_vector (GfsVariable ** v, guint n)
+{
+ guint i, j;
+
+ g_return_if_fail (v != NULL);
+ g_return_if_fail (n <= FTT_DIMENSION);
+
+ for (i = 0; i < n; i++) {
+ g_return_if_fail (v[i] != NULL);
+ v[i]->component = i;
+ for (j = 0; j < n; j++)
+ v[i]->vector[j] = v[j];
+ }
+}
+
/* GfsVariableTracer: object */
static void variable_tracer_read (GtsObject ** o, GtsFile * fp)
diff --git a/src/variable.h b/src/variable.h
index 43d60c7..fa57df3 100644
--- a/src/variable.h
+++ b/src/variable.h
@@ -40,6 +40,7 @@ struct _GfsVariable {
/*< public >*/
guint i;
FttComponent component;
+ GfsVariable * vector[FTT_DIMENSION];
gchar * name, * description;
gboolean centered;
GfsVariableFineCoarseFunc fine_coarse, coarse_fine;
@@ -84,7 +85,8 @@ GSList * gfs_variables_from_list (GSList * i,
gchar ** error);
void gfs_variables_swap (GfsVariable * v1,
GfsVariable * v2);
-#define gfs_variable_set_vector(v, c) ((v)->component = (c))
+void gfs_variable_set_vector (GfsVariable ** v,
+ guint n);
/* GfsVariableTracer: header */
--
Gerris Flow Solver
More information about the debian-science-commits
mailing list