Skip to content

Commit a2ce6a8

Browse files
committed
Add grad-div Vector FE Ex. 10
1 parent 2a4aedf commit a2ce6a8

9 files changed

Lines changed: 729 additions & 2 deletions

File tree

configure.ac

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -480,6 +480,7 @@ AS_IF([test "x$enableexamples" = "xyes"],
480480
examples/vector_fe/vector_fe_ex7/Makefile
481481
examples/vector_fe/vector_fe_ex8/Makefile
482482
examples/vector_fe/vector_fe_ex9/Makefile
483+
examples/vector_fe/vector_fe_ex10/Makefile
483484
examples/Makefile])
484485
])
485486

examples/Makefile.am

Lines changed: 2 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -74,7 +74,8 @@ SUBDIRS = \
7474
vector_fe/vector_fe_ex6 \
7575
vector_fe/vector_fe_ex7 \
7676
vector_fe/vector_fe_ex8 \
77-
vector_fe/vector_fe_ex9
77+
vector_fe/vector_fe_ex9 \
78+
vector_fe/vector_fe_ex10
7879

7980
AUTOMAKE_OPTIONS = subdir-objects
8081

Lines changed: 21 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,21 @@
1+
example_name = vector_fe_ex10
2+
check_SCRIPTS = run.sh
3+
install_dir = $(examples_install_path)/vector_fe/ex10
4+
data = grad_div_exact_solution.h solution_function.h vector_fe_ex10.C vector_fe_ex10.in run.sh
5+
sources = $(data)
6+
7+
CLEANFILES = out.e
8+
9+
# also need to link files for VPATH builds
10+
if LIBMESH_VPATH_BUILD
11+
BUILT_SOURCES = .linkstamp
12+
.linkstamp:
13+
-rm -f vector_fe_ex10.in && $(LN_S) -f $(srcdir)/vector_fe_ex10.in .
14+
$(AM_V_GEN)touch .linkstamp
15+
16+
CLEANFILES += vector_fe_ex10.in .linkstamp
17+
endif
18+
19+
##############################################
20+
# include common example environment
21+
include $(top_srcdir)/examples/Make.common
Lines changed: 63 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,63 @@
1+
// The libMesh Finite Element Library.
2+
// Copyright (C) 2002-2024 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner
3+
4+
// This library is free software; you can redistribute it and/or
5+
// modify it under the terms of the GNU Lesser General Public
6+
// License as published by the Free Software Foundation; either
7+
// version 2.1 of the License, or (at your option) any later version.
8+
9+
// This library is distributed in the hope that it will be useful,
10+
// but WITHOUT ANY WARRANTY; without even the implied warranty of
11+
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
12+
// Lesser General Public License for more details.
13+
14+
// You should have received a copy of the GNU Lesser General Public
15+
// License along with this library; if not, write to the Free Software
16+
// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
17+
18+
#ifndef GRAD_DIV_EXACT_SOLUTION_H
19+
#define GRAD_DIV_EXACT_SOLUTION_H
20+
21+
#include "libmesh/libmesh_common.h"
22+
#include "libmesh/vector_value.h"
23+
24+
using namespace libMesh;
25+
26+
class GradDivExactSolution
27+
{
28+
public:
29+
GradDivExactSolution() = default;
30+
~GradDivExactSolution() = default;
31+
32+
RealGradient operator() (Real x, Real y, Real z)
33+
{
34+
libmesh_ignore(z);
35+
36+
const Real ux = cos(k*x)*sin(k*y);
37+
const Real uy = sin(k*x)*cos(k*y);
38+
39+
return RealGradient(ux, uy);
40+
}
41+
42+
RealTensor grad(Real x, Real y, Real z)
43+
{
44+
libmesh_ignore(z);
45+
46+
const Real dux_dx = -k*sin(k*x)*sin(k*y);
47+
const Real dux_dy = k*cos(k*x)*cos(k*y);
48+
const Real duy_dx = dux_dy;
49+
const Real duy_dy = dux_dx;
50+
51+
return RealTensor(dux_dx, dux_dy, Real(0), duy_dx, duy_dy);
52+
}
53+
54+
RealGradient forcing(Real x, Real y, Real z)
55+
{
56+
return (2*k*k + 1)*operator()(x, y, z);
57+
}
58+
59+
private:
60+
const Real k = pi;
61+
};
62+
63+
#endif // GRAD_DIV_EXACT_SOLUTION_H
Lines changed: 41 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,41 @@
1+
#!/bin/sh
2+
3+
#set -x
4+
5+
. "$LIBMESH_DIR"/examples/run_common.sh
6+
7+
example_name=vector_fe_ex10
8+
9+
# Note: these problems are particularly ill-conditioned, so we currently
10+
# resort to a direct solver.
11+
12+
options="dim=2 element_type=TRI6 -pc_type lu"
13+
run_example "$example_name" "$options"
14+
15+
options="dim=2 element_type=TRI7 -pc_type lu"
16+
run_example "$example_name" "$options"
17+
18+
options="dim=2 element_type=QUAD8 -pc_type lu"
19+
run_example "$example_name" "$options"
20+
21+
options="dim=2 element_type=QUAD9 -pc_type lu"
22+
run_example "$example_name" "$options"
23+
24+
options="dim=2 order=3 element_type=TRI6 -pc_type lu"
25+
run_example "$example_name" "$options"
26+
27+
options="dim=2 order=3 element_type=TRI7 -pc_type lu"
28+
run_example "$example_name" "$options"
29+
30+
options="dim=2 order=3 element_type=QUAD8 -pc_type lu"
31+
run_example "$example_name" "$options"
32+
33+
options="dim=2 order=3 element_type=QUAD9 -pc_type lu"
34+
run_example "$example_name" "$options"
35+
36+
# Subdividing each hex into 24 tets gets expensive in dbg...
37+
options="dim=3 element_type=TET14 grid_size=6 -pc_type lu"
38+
run_example "$example_name" "$options"
39+
40+
options="dim=3 element_type=HEX27 -pc_type lu"
41+
run_example "$example_name" "$options"
Lines changed: 109 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,109 @@
1+
// The libMesh Finite Element Library.
2+
// Copyright (C) 2002-2024 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner
3+
4+
// This library is free software; you can redistribute it and/or
5+
// modify it under the terms of the GNU Lesser General Public
6+
// License as published by the Free Software Foundation; either
7+
// version 2.1 of the License, or (at your option) any later version.
8+
9+
// This library is distributed in the hope that it will be useful,
10+
// but WITHOUT ANY WARRANTY; without even the implied warranty of
11+
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
12+
// Lesser General Public License for more details.
13+
14+
// You should have received a copy of the GNU Lesser General Public
15+
// License along with this library; if not, write to the Free Software
16+
// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
17+
18+
#ifndef SOLUTION_FUNCTION_H
19+
#define SOLUTION_FUNCTION_H
20+
21+
// libMesh includes
22+
#include "libmesh/function_base.h"
23+
24+
// Example includes
25+
#include "grad_div_exact_solution.h"
26+
27+
// C++ includes
28+
#include <memory>
29+
30+
using namespace libMesh;
31+
32+
class SolutionFunction : public FunctionBase<Number>
33+
{
34+
public:
35+
36+
SolutionFunction() = default;
37+
~SolutionFunction() = default;
38+
39+
virtual Number operator() (const Point &,
40+
const Real = 0)
41+
{ libmesh_not_implemented(); }
42+
43+
virtual void operator() (const Point & p,
44+
const Real,
45+
DenseVector<Number> & output)
46+
{
47+
output.zero();
48+
const Real x=p(0), y=p(1), z=p(2);
49+
output(0) = soln(x, y, z)(0);
50+
output(1) = soln(x, y, z)(1);
51+
output(2) = soln(x, y, z)(2);
52+
}
53+
54+
virtual Number component(unsigned int component_in,
55+
const Point & p,
56+
const Real)
57+
{
58+
DenseVector<Number> outvec(3);
59+
(*this)(p, 0, outvec);
60+
return outvec(component_in);
61+
}
62+
63+
virtual std::unique_ptr<FunctionBase<Number>> clone() const
64+
{ return std::make_unique<SolutionFunction>(); }
65+
66+
private:
67+
68+
GradDivExactSolution soln;
69+
};
70+
71+
class SolutionGradient : public FunctionBase<Gradient>
72+
{
73+
public:
74+
75+
SolutionGradient() = default;
76+
~SolutionGradient() = default;
77+
78+
virtual Gradient operator() (const Point &, const Real = 0)
79+
{ libmesh_not_implemented(); }
80+
81+
virtual void operator() (const Point & p,
82+
const Real,
83+
DenseVector<Gradient> & output)
84+
{
85+
output.zero();
86+
const Real x=p(0), y=p(1), z=p(2);
87+
output(0) = soln.grad(x, y, z).row(0);
88+
output(1) = soln.grad(x, y, z).row(1);
89+
output(2) = soln.grad(x, y, z).row(2);
90+
}
91+
92+
virtual Gradient component(unsigned int component_in,
93+
const Point & p,
94+
const Real)
95+
{
96+
DenseVector<Gradient> outvec(3);
97+
(*this)(p, 0, outvec);
98+
return outvec(component_in);
99+
}
100+
101+
virtual std::unique_ptr<FunctionBase<Gradient>> clone() const
102+
{ return std::make_unique<SolutionGradient>(); }
103+
104+
private:
105+
106+
GradDivExactSolution soln;
107+
};
108+
109+
#endif // SOLUTION_FUNCTION_H

0 commit comments

Comments
 (0)