Skip to content

Commit 5f7cb5f

Browse files
committed
ExternalMeshDeformation class
1 parent 29c1e20 commit 5f7cb5f

File tree

10 files changed

+924
-14
lines changed

10 files changed

+924
-14
lines changed

include/aspect/mesh_deformation/external_tool_interface.h

Lines changed: 6 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -112,6 +112,11 @@ namespace aspect
112112
AffineConstraints<double> &mesh_velocity_constraints,
113113
const std::set<types::boundary_id> &boundary_ids) const override;
114114

115+
/**
116+
* Given all solution variables at each surface point, compute velocities at these points.
117+
*
118+
* This needs to be implemented by the derived class and implement the "surface evolution".
119+
*/
115120
virtual
116121
std::vector<Tensor<1,dim>>
117122
compute_updated_velocities_at_points (const std::vector<std::vector<double>> &current_solution_at_points) const = 0;
@@ -176,8 +181,7 @@ namespace aspect
176181
LinearAlgebra::Vector
177182
interpolate_velocities_to_surface_points (const std::vector<Tensor<1,dim>> &vertical_velocities) const;
178183

179-
180-
private:
184+
protected:
181185

182186
std::vector<Point<dim>> evaluation_points;
183187

source/mesh_deformation/external_tool_interface.cc

Lines changed: 380 additions & 12 deletions
Large diffs are not rendered by default.
Lines changed: 136 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,136 @@
1+
/*
2+
Copyright (C) 2025 - 2025 by the authors of the ASPECT code.
3+
4+
This file is part of ASPECT.
5+
6+
ASPECT is free software; you can redistribute it and/or modify
7+
it under the terms of the GNU General Public License as published by
8+
the Free Software Foundation; either version 2, or (at your option)
9+
any later version.
10+
11+
ASPECT is distributed in the hope that it will be useful,
12+
but WITHOUT ANY WARRANTY; without even the implied warranty of
13+
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14+
GNU General Public License for more details.
15+
16+
You should have received a copy of the GNU General Public License
17+
along with ASPECT; see the file LICENSE. If not see
18+
<http://www.gnu.org/licenses/>.
19+
*/
20+
21+
#include <aspect/simulator_signals.h>
22+
#include <aspect/simulator_access.h>
23+
24+
#include <iostream>
25+
26+
using namespace aspect;
27+
28+
#include <aspect/mesh_deformation/interface.h>
29+
#include <aspect/gravity_model/interface.h>
30+
#include <aspect/simulator_access.h>
31+
#include <aspect/mesh_deformation/external_tool_interface.h>
32+
33+
34+
namespace aspect
35+
{
36+
namespace MeshDeformation
37+
{
38+
template <int dim>
39+
class TestExternalDeformation : public ExternalToolInterface<dim>
40+
{
41+
public:
42+
/**
43+
* Constructor.
44+
*/
45+
TestExternalDeformation() = default;
46+
47+
void initialize() override
48+
{
49+
this->get_pcout() << "initialize()" << std::endl;
50+
}
51+
52+
void
53+
update () override
54+
{
55+
this->get_pcout() << "update()" << std::endl;
56+
57+
if (!this->remote_point_evaluator)
58+
{
59+
std::vector<Point<dim>> points;
60+
this->get_pcout() << "\tsetting points" << std::endl;
61+
62+
if (Utilities::MPI::this_mpi_process(this->get_mpi_communicator())==0)
63+
{
64+
if constexpr (dim == 2)
65+
{
66+
points.emplace_back(0.3, 1.0);
67+
points.emplace_back(0.5, 1.0);
68+
}
69+
else
70+
{
71+
points.emplace_back(0.3, 0.7, 1.0);
72+
points.emplace_back(0.5, 0.5, 1.0);
73+
points.emplace_back(0.7, 0.2, 1.0);
74+
points.emplace_back(0.7, 0.91, 1.0);
75+
}
76+
}
77+
else
78+
{
79+
if constexpr (dim == 2)
80+
{
81+
points.emplace_back(0.72, 1.0);
82+
}
83+
else
84+
{
85+
points.emplace_back(0.9, 0.9, 1.0);
86+
}
87+
}
88+
this->set_evaluation_points (points);
89+
}
90+
}
91+
92+
virtual
93+
std::vector<Tensor<1,dim>>
94+
compute_updated_velocities_at_points (const std::vector<std::vector<double>> &current_solution_at_points) const override
95+
{
96+
this->get_pcout() << "compute_updated_velocities_at_points()" << std::endl;
97+
98+
Assert(current_solution_at_points.size() == this->evaluation_points.size(), ExcInternalError());
99+
std::vector<Tensor<1,dim>> velocities(current_solution_at_points.size(), Tensor<1,dim>());
100+
if (velocities.size()>1)
101+
{
102+
velocities[0][dim-1]=30.0;
103+
velocities[1][dim-1]=-3.5;
104+
}
105+
return velocities;
106+
}
107+
108+
109+
110+
Tensor<1,dim>
111+
compute_initial_deformation_on_boundary(const types::boundary_id /*boundary_indicator*/,
112+
const Point<dim> &position) const override
113+
{
114+
const Tensor<1,dim> gravity = this->get_gravity_model().gravity_vector(position);
115+
Tensor<1,dim> topography_direction;
116+
if (gravity.norm() > 0.0)
117+
topography_direction = -gravity / gravity.norm();
118+
119+
const double topography_amplitude = (position[0]>=0.5) ? (0.05 * (1.+std::cos(2.*numbers::PI*position[0]))) : 0.0;
120+
return topography_amplitude * topography_direction;
121+
}
122+
};
123+
}
124+
}
125+
126+
127+
// explicit instantiation of the functions we implement in this file
128+
namespace aspect
129+
{
130+
namespace MeshDeformation
131+
{
132+
ASPECT_REGISTER_MESH_DEFORMATION_MODEL(TestExternalDeformation,
133+
"external deformation",
134+
"")
135+
}
136+
}
Lines changed: 96 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,96 @@
1+
# Test the ExternalToolInterface class for mesh deformation
2+
# derived class written as a plugin, 2d box geometry, initial
3+
# topography coming from the plugin (and not initial topography)
4+
#
5+
# MPI: 2
6+
7+
set Dimension = 2
8+
set Use years in output instead of seconds = false
9+
set End time = 0.004
10+
set Maximum time step = 0.0005
11+
set Nonlinear solver scheme = no Advection, no Stokes
12+
set Pressure normalization = surface
13+
set Surface pressure = 0
14+
15+
subsection Geometry model
16+
set Model name = box
17+
18+
subsection Box
19+
set X extent = 1
20+
set Y extent = 1
21+
end
22+
end
23+
24+
# Temperature effects are ignored
25+
subsection Initial temperature model
26+
set Model name = function
27+
28+
subsection Function
29+
set Function expression = 0
30+
end
31+
end
32+
33+
subsection Boundary temperature model
34+
set Fixed temperature boundary indicators = bottom, top
35+
set List of model names = initial temperature
36+
end
37+
38+
# Free slip on all boundaries
39+
subsection Boundary velocity model
40+
set Tangential velocity boundary indicators = left, right, bottom, top
41+
end
42+
43+
subsection Mesh deformation
44+
set Mesh deformation boundary indicators = top: external deformation
45+
set Additional tangential mesh velocity boundary indicators = left, right
46+
end
47+
48+
# Vertical gravity
49+
subsection Gravity model
50+
set Model name = vertical
51+
52+
subsection Vertical
53+
set Magnitude = 1
54+
end
55+
end
56+
57+
# One material with unity properties
58+
subsection Material model
59+
set Model name = simple
60+
61+
subsection Simple model
62+
set Reference density = 1
63+
set Reference specific heat = 1
64+
set Reference temperature = 0
65+
set Thermal conductivity = 1
66+
set Thermal expansion coefficient = 1
67+
set Viscosity = 1
68+
end
69+
end
70+
71+
# We also have to specify that we want to use the Boussinesq
72+
# approximation (assuming the density in the temperature
73+
# equation to be constant, and incompressibility).
74+
subsection Formulation
75+
set Formulation = Boussinesq approximation
76+
end
77+
78+
# We here use a globally refined mesh without
79+
# adaptive mesh refinement.
80+
subsection Mesh refinement
81+
set Initial global refinement = 4
82+
set Initial adaptive refinement = 0
83+
set Time steps between mesh refinement = 0
84+
end
85+
86+
subsection Postprocess
87+
set List of postprocessors = visualization, topography
88+
89+
subsection Visualization
90+
set Time between graphical output = 0
91+
set Output mesh velocity = true
92+
set Output mesh displacement = true
93+
set Output undeformed mesh = false
94+
set Interpolate output = false
95+
end
96+
end
Lines changed: 105 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,105 @@
1+
-----------------------------------------------------------------------------
2+
-----------------------------------------------------------------------------
3+
-----------------------------------------------------------------------------
4+
5+
Loading shared library <./libmesh_deformation_external_01.debug.so>
6+
7+
initialize()
8+
-----------------------------------------------------------------------------
9+
-----------------------------------------------------------------------------
10+
Number of active cells: 256 (on 5 levels)
11+
Number of degrees of freedom: 3,556 (2,178+289+1,089)
12+
13+
Number of mesh deformation degrees of freedom: 578
14+
Solving mesh displacement system... 6 iterations.
15+
*** Timestep 0: t=0 seconds, dt=0 seconds
16+
update()
17+
setting points
18+
compute_updated_velocities_at_points()
19+
Solving mesh displacement system... 4 iterations.
20+
21+
Postprocessing:
22+
Writing graphical output: output-mesh_deformation_external_01/solution/solution-00000
23+
Topography min/max: 0 m, 0.1 m
24+
25+
*** Timestep 1: t=0.0005 seconds, dt=0.0005 seconds
26+
update()
27+
compute_updated_velocities_at_points()
28+
Solving mesh displacement system... 4 iterations.
29+
30+
Postprocessing:
31+
Writing graphical output: output-mesh_deformation_external_01/solution/solution-00001
32+
Topography min/max: -0.00175 m, 0.1 m
33+
34+
*** Timestep 2: t=0.001 seconds, dt=0.0005 seconds
35+
update()
36+
compute_updated_velocities_at_points()
37+
Solving mesh displacement system... 4 iterations.
38+
39+
Postprocessing:
40+
Writing graphical output: output-mesh_deformation_external_01/solution/solution-00002
41+
Topography min/max: -0.0035 m, 0.1 m
42+
43+
*** Timestep 3: t=0.0015 seconds, dt=0.0005 seconds
44+
update()
45+
compute_updated_velocities_at_points()
46+
Solving mesh displacement system... 4 iterations.
47+
48+
Postprocessing:
49+
Writing graphical output: output-mesh_deformation_external_01/solution/solution-00003
50+
Topography min/max: -0.00525 m, 0.1 m
51+
52+
*** Timestep 4: t=0.002 seconds, dt=0.0005 seconds
53+
update()
54+
compute_updated_velocities_at_points()
55+
Solving mesh displacement system... 4 iterations.
56+
57+
Postprocessing:
58+
Writing graphical output: output-mesh_deformation_external_01/solution/solution-00004
59+
Topography min/max: -0.007 m, 0.1 m
60+
61+
*** Timestep 5: t=0.0025 seconds, dt=0.0005 seconds
62+
update()
63+
compute_updated_velocities_at_points()
64+
Solving mesh displacement system... 4 iterations.
65+
66+
Postprocessing:
67+
Writing graphical output: output-mesh_deformation_external_01/solution/solution-00005
68+
Topography min/max: -0.00875 m, 0.1 m
69+
70+
*** Timestep 6: t=0.003 seconds, dt=0.0005 seconds
71+
update()
72+
compute_updated_velocities_at_points()
73+
Solving mesh displacement system... 4 iterations.
74+
75+
Postprocessing:
76+
Writing graphical output: output-mesh_deformation_external_01/solution/solution-00006
77+
Topography min/max: -0.0105 m, 0.1 m
78+
79+
*** Timestep 7: t=0.0035 seconds, dt=0.0005 seconds
80+
update()
81+
compute_updated_velocities_at_points()
82+
Solving mesh displacement system... 4 iterations.
83+
84+
Postprocessing:
85+
Writing graphical output: output-mesh_deformation_external_01/solution/solution-00007
86+
Topography min/max: -0.01225 m, 0.105 m
87+
88+
*** Timestep 8: t=0.004 seconds, dt=0.0005 seconds
89+
update()
90+
compute_updated_velocities_at_points()
91+
Solving mesh displacement system... 4 iterations.
92+
93+
Postprocessing:
94+
Writing graphical output: output-mesh_deformation_external_01/solution/solution-00008
95+
Topography min/max: -0.014 m, 0.12 m
96+
97+
Termination requested by criterion: end time
98+
99+
100+
+---------------------------------------------+------------+------------+
101+
+---------------------------------+-----------+------------+------------+
102+
+---------------------------------+-----------+------------+------------+
103+
104+
-----------------------------------------------------------------------------
105+
-----------------------------------------------------------------------------
Lines changed: 19 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,19 @@
1+
# 1: Time step number
2+
# 2: Time (seconds)
3+
# 3: Time step size (seconds)
4+
# 4: Number of mesh cells
5+
# 5: Number of Stokes degrees of freedom
6+
# 6: Number of temperature degrees of freedom
7+
# 7: Number of nonlinear iterations
8+
# 8: Visualization file name
9+
# 9: Minimum topography (m)
10+
# 10: Maximum topography (m)
11+
0 0.000000000000e+00 0.000000000000e+00 256 2467 1089 0 output-mesh_deformation_external_01/solution/solution-00000 0.00000000e+00 1.00000000e-01
12+
1 5.000000000000e-04 5.000000000000e-04 256 2467 1089 0 output-mesh_deformation_external_01/solution/solution-00001 -1.75000000e-03 1.00000000e-01
13+
2 1.000000000000e-03 5.000000000000e-04 256 2467 1089 0 output-mesh_deformation_external_01/solution/solution-00002 -3.50000000e-03 1.00000000e-01
14+
3 1.500000000000e-03 5.000000000000e-04 256 2467 1089 0 output-mesh_deformation_external_01/solution/solution-00003 -5.25000000e-03 1.00000000e-01
15+
4 2.000000000000e-03 5.000000000000e-04 256 2467 1089 0 output-mesh_deformation_external_01/solution/solution-00004 -7.00000000e-03 1.00000000e-01
16+
5 2.500000000000e-03 5.000000000000e-04 256 2467 1089 0 output-mesh_deformation_external_01/solution/solution-00005 -8.75000000e-03 1.00000000e-01
17+
6 3.000000000000e-03 5.000000000000e-04 256 2467 1089 0 output-mesh_deformation_external_01/solution/solution-00006 -1.05000000e-02 1.00000000e-01
18+
7 3.500000000000e-03 5.000000000000e-04 256 2467 1089 0 output-mesh_deformation_external_01/solution/solution-00007 -1.22500000e-02 1.05000000e-01
19+
8 4.000000000000e-03 5.000000000000e-04 256 2467 1089 0 output-mesh_deformation_external_01/solution/solution-00008 -1.40000000e-02 1.20000000e-01

0 commit comments

Comments
 (0)