[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
[Getfem-commits] (no subject)
From: |
Tetsuo Koyama |
Subject: |
[Getfem-commits] (no subject) |
Date: |
Sat, 1 Jun 2019 06:22:41 -0400 (EDT) |
branch: devel-tetsuo-houbolt
commit 398315140ef3cb2913827af1c4bb3c784539efca
Author: Tetsuo Koyama <address@hidden>
Date: Sat Jun 1 19:20:19 2019 +0900
Add Houbolt method
---
src/getfem/getfem_models.h | 2 +
src/getfem_models.cc | 112 +++++++++++++++++++++++++++++++++++++++++++++
2 files changed, 114 insertions(+)
diff --git a/src/getfem/getfem_models.h b/src/getfem/getfem_models.h
index 697b5fa..9d10519 100644
--- a/src/getfem/getfem_models.h
+++ b/src/getfem/getfem_models.h
@@ -1180,6 +1180,8 @@ namespace getfem {
void add_Newmark_scheme(model &md, const std::string &varname,
scalar_type beta, scalar_type gamma);
+ void add_Houbolt_scheme(model &md, const std::string &varname);
+
//=========================================================================
diff --git a/src/getfem_models.cc b/src/getfem_models.cc
index d6b7967..87c9e80 100644
--- a/src/getfem_models.cc
+++ b/src/getfem_models.cc
@@ -1537,6 +1537,118 @@ namespace getfem {
md.add_time_integration_scheme(varname, ptsc);
}
+ // ----------------------------------------------------------------------
+ //
+ // Houbolt method
+ //
+ // ----------------------------------------------------------------------
+
+ class APIDECL Houbolt_scheme
+ : public virtual_time_scheme {
+
+ std::string U, U01, U02, U03, V, A;
+
+ public:
+ // V = 1/(6*dt)*(11*U-18*U01+9*U02-2*U03)
+ // A = 1/(dt**2)*(2*U-5*U01+4*U02-U03)
+ virtual void init_affine_dependent_variables(model &md) const {
+ scalar_type dt = md.get_time_step();
+ scalar_type a0 = scalar_type(2)/(dt*dt);
+ scalar_type a1 = scalar_type(5)/(dt*dt);
+ scalar_type a2 = scalar_type(4)/(dt*dt);
+ scalar_type a3 = scalar_type(1)/(dt*dt);
+ scalar_type b0 = scalar_type(11)/(scalar_type(6)*dt);
+ scalar_type b1 = scalar_type(18)/(scalar_type(6)*dt);
+ scalar_type b2 = scalar_type(9)/(scalar_type(6)*dt);
+ scalar_type b3 = scalar_type(2)/(scalar_type(6)*dt);
+
+ md.set_factor_of_variable(V, b0);
+ md.set_factor_of_variable(A, a0);
+ if (md.is_complex()) {
+ gmm::add(gmm::scaled(md.complex_variable(U01), -complex_type(b1)),
+ gmm::scaled(md.complex_variable(U02), complex_type(b2)),
+ md.set_complex_constant_part(V));
+ gmm::add(gmm::scaled(md.complex_variable(U03), -complex_type(b3)),
+ md.set_complex_constant_part(V));
+ gmm::add(gmm::scaled(md.complex_variable(U01), -complex_type(a1)),
+ gmm::scaled(md.complex_variable(U02), complex_type(a2)),
+ md.set_complex_constant_part(A));
+ gmm::add(gmm::scaled(md.complex_variable(U03), -complex_type(a3)),
+ md.set_complex_constant_part(A));
+ } else {
+ gmm::add(gmm::scaled(md.real_variable(U01), -b1),
+ gmm::scaled(md.real_variable(U02), b2),
+ md.set_real_constant_part(V));
+ gmm::add(gmm::scaled(md.real_variable(U03), -b3),
+ md.set_real_constant_part(V));
+ gmm::add(gmm::scaled(md.real_variable(U01), -a1),
+ gmm::scaled(md.real_variable(U02), a2),
+ md.set_real_constant_part(A));
+ gmm::add(gmm::scaled(md.real_variable(U03), -a3),
+ md.set_real_constant_part(A));
+ }
+ }
+
+ virtual void init_affine_dependent_variables_precomputation(model &md)
+ const {
+ (void) md;
+ }
+
+ virtual void time_derivative_to_be_initialized
+ (std::string &name_v, std::string &name_previous_v) const {
+ (void) name_v;
+ (void) name_previous_v;
+ }
+
+ virtual void shift_variables(model &md) const {
+ if (md.is_complex()) {
+ gmm::copy(md.complex_variable(U02), md.set_complex_variable(U03));
+ gmm::copy(md.complex_variable(U01), md.set_complex_variable(U02));
+ gmm::copy(md.complex_variable(U), md.set_complex_variable(U01));
+ } else {
+ gmm::copy(md.real_variable(U02), md.set_complex_variable(U03));
+ gmm::copy(md.real_variable(U01), md.set_complex_variable(U02));
+ gmm::copy(md.real_variable(U), md.set_complex_variable(U01));
+ }
+ }
+
+
+ Houbolt_scheme(model &md, std::string varname) {
+ U = varname;
+ U01 = "Previous_" + U;
+ U02 = "Previous2_" + U;
+ U03 = "Previous3_" + U;
+ V = "Dot_" + U;
+ A = "Dot2_" + U;
+
+ if (!(md.variable_exists(V)))
+ md.add_affine_dependent_variable(V, U);
+ if (!(md.variable_exists(A)))
+ md.add_affine_dependent_variable(A, U);
+
+ const mesh_fem *mf = md.pmesh_fem_of_variable(U);
+ size_type s = md.is_complex() ? gmm::vect_size(md.complex_variable(U))
+ : gmm::vect_size(md.real_variable(U));
+
+ if (mf) {
+ if (!(md.variable_exists(U01))) md.add_fem_data(U01, *mf);
+ if (!(md.variable_exists(U02))) md.add_fem_data(U02, *mf);
+ if (!(md.variable_exists(U03))) md.add_fem_data(U03, *mf);
+ } else {
+ if (!(md.variable_exists(U01))) md.add_fixed_size_data(U01, s);
+ if (!(md.variable_exists(U02))) md.add_fixed_size_data(U02, s);
+ if (!(md.variable_exists(U03))) md.add_fixed_size_data(U03, s);
+ }
+
+ }
+
+ };
+
+ void add_Houbolt_scheme(model &md, const std::string &varname) {
+ ptime_scheme ptsc = std::make_shared<Houbolt_scheme>
+ (md, varname);
+ md.add_time_integration_scheme(varname, ptsc);
+ }