From 905dafa9a559cb33af055d90a37e2c3767ae2254 Mon Sep 17 00:00:00 2001
From: Elias Pipping <elias.pipping@fu-berlin.de>
Date: Tue, 11 Sep 2012 15:51:31 +0200
Subject: [PATCH] Solve velocity problem w/ corresponding Dirichlet
---
src/one-body-sample.py | 2 +-
src/timestepping.cc | 37 +++++++++++++------------------------
2 files changed, 14 insertions(+), 25 deletions(-)
diff --git a/src/one-body-sample.py b/src/one-body-sample.py
index 77b2f491..be879d3c 100644
--- a/src/one-body-sample.py
+++ b/src/one-body-sample.py
@@ -13,7 +13,7 @@ class neumannCondition:
class dirichletCondition:
def __call__(self, x):
- return x * .005
+ return .005
Functions = {
'neumannCondition' : neumannCondition(),
diff --git a/src/timestepping.cc b/src/timestepping.cc
index 98b89e79..dd9ad0e0 100644
--- a/src/timestepping.cc
+++ b/src/timestepping.cc
@@ -56,23 +56,20 @@ class ImplicitEuler
problem_A = this->A;
problem_A *= this->tau;
- // Treat as Delta u_n for now
- problem_iterate = this->u_old;
- if (u_old_old_ptr)
+ // Use the old velocity as an initial iterate if possible
+ if (u_old_old_ptr) {
+ problem_iterate = this->u_old;
problem_iterate -= *u_old_old_ptr;
+ problem_iterate /= this->tau;
+ } else
+ problem_iterate = 0.0;
for (size_t i = 0; i < this->dirichletNodes.size(); ++i)
if (this->dirichletNodes[i].count() == dim) {
- double val;
- this->dirichletFunction.evaluate(this->time, val);
- problem_iterate[i] = 0; // Everything prescribed
- problem_iterate[i][0] =
- val - this->u_old[i][0]; // Time-dependent X direction
+ problem_iterate[i] = 0;
+ this->dirichletFunction.evaluate(this->time, problem_iterate[i][0]);
} else if (this->dirichletNodes[i][1])
- problem_iterate[i][1] = 0; // Y direction described
-
- // Make it a velocity
- problem_iterate /= this->tau;
+ problem_iterate[i][1] = 0; // Y direction prescribed
}
void virtual extractSolution(VectorType const &problem_iterate,
@@ -112,25 +109,17 @@ class ImplicitTwoStep
problem_A = this->A;
problem_A *= 2.0 / 3.0 * this->tau;
- // The finite difference makes a good start (treat it as such for now)
+ // Use the old velocity as an initial iterate
problem_iterate = this->u_old;
problem_iterate -= *u_old_old_ptr;
+ problem_iterate /= this->tau;
for (size_t i = 0; i < this->dirichletNodes.size(); ++i)
if (this->dirichletNodes[i].count() == dim) {
- double val;
- this->dirichletFunction.evaluate(this->time, val);
problem_iterate[i] = 0;
- problem_iterate[i].axpy(-2, this->u_old[i]);
- problem_iterate[i].axpy(.5, (*u_old_old_ptr)[i]);
- problem_iterate[i][0] =
- 1.5 * val - 2 * this->u_old[i][0] + .5 * (*u_old_old_ptr)[i][0];
+ this->dirichletFunction.evaluate(this->time, problem_iterate[i][0]);
} else if (this->dirichletNodes[i][1])
- // Y direction described
- problem_iterate[i][1] =
- -2 * this->u_old[i][1] + .5 * (*u_old_old_ptr)[i][1];
- // Now treat it as a velocity
- problem_iterate /= this->tau;
+ problem_iterate[i][1] = 0; // Y direction prescribed
}
void virtual extractSolution(VectorType const &problem_iterate,
--
GitLab