[cig-commits] r13857 - in cs/cigma/trunk/examples: . laplace1 laplace2
luis at geodynamics.org
luis at geodynamics.org
Tue Jan 13 03:11:59 PST 2009
Author: luis
Date: 2009-01-13 03:11:59 -0800 (Tue, 13 Jan 2009)
New Revision: 13857
Added:
cs/cigma/trunk/examples/laplace1/
cs/cigma/trunk/examples/laplace1/Makefile
cs/cigma/trunk/examples/laplace1/README
cs/cigma/trunk/examples/laplace1/laplace.cc
cs/cigma/trunk/examples/laplace2/
cs/cigma/trunk/examples/laplace2/Makefile
cs/cigma/trunk/examples/laplace2/README
cs/cigma/trunk/examples/laplace2/exact_solution.C
cs/cigma/trunk/examples/laplace2/exact_solution.H
cs/cigma/trunk/examples/laplace2/laplace.C
Modified:
cs/cigma/trunk/examples/README
Log:
Added source code for a simple example. We output the solutions
in a format easily understood by cigma.
Specifically, we use a simple Poisson problem, solving in two different
ways and using different elements, to make sure we can use cigma to
calculate the right order of convergence of each solution method.
The first example uses the deal.II library, while the second
uses the libMesh library.
Modified: cs/cigma/trunk/examples/README
===================================================================
--- cs/cigma/trunk/examples/README 2009-01-13 00:03:43 UTC (rev 13856)
+++ cs/cigma/trunk/examples/README 2009-01-13 11:11:59 UTC (rev 13857)
@@ -1 +1,2 @@
-This directory contains examples
+In this directory you will find source code for reproducing most
+of the examples discussed in the Cigma manual.
Added: cs/cigma/trunk/examples/laplace1/Makefile
===================================================================
--- cs/cigma/trunk/examples/laplace1/Makefile (rev 0)
+++ cs/cigma/trunk/examples/laplace1/Makefile 2009-01-13 11:11:59 UTC (rev 13857)
@@ -0,0 +1,64 @@
+###############################################################################
+# -*- Makefile -*-
+#
+
+D = $(HOME)/dealii
+
+target = laplace
+target-source = laplace.cc
+
+###############################################################################
+
+debug-mode = on
+clean-up-files = *gmv *gnuplot *gpl *eps *pov
+target.exe = $(target)$(EXEEXT)
+
+include $D/common/Make.global_options
+
+libs.g = \
+ $(lib-deal2-1d.g) \
+ $(lib-deal2-2d.g) \
+ $(lib-deal2-3d.g) \
+ $(lib-lac.g) \
+ $(lib-base.g)
+
+libs.o = \
+ $(lib-deal2-1d.o) \
+ $(lib-deal2-2d.o) \
+ $(lib-deal2-3d.o) \
+ $(lib-lac.o) \
+ $(lib-base.o)
+
+ifeq ($(debug-mode),on)
+ libraries = $(target).g.$(OBJEXT) $(libs.g)
+else
+ libraries = $(target).$(OBJEXT) $(libs.o)
+endif
+
+$(target.exe): $(libraries)
+ $(CXX) -o $@ $^ $(LIBS) $(LDFLAGS)
+
+run: $(target.exe)
+ @echo === running $<
+ @./$(target.exe)
+
+clean:
+ -rm -f *.$(OBJEXT) *~ Makefile.dep $(target.exe) $(clean-up-files)
+
+./%.g.$(OBJEXT):
+ @echo === debug $(<F)
+ $(CXX) $(CXXFLAGS.g) -c $< -o $@
+
+./%.$(OBJEXT):
+ @echo === optimized $(<F)
+ $(CXX) $(CXXFLAGS.o) -c $< -o $@
+
+.PHONY: run clean
+
+Makefile.dep: $(target-source) Makefile $(shell echo $D/*/include/*/*.h)
+ @echo === remaking $@
+ @$D/common/scripts/make_dependencies $(INCLUDE) -B. $(target-source) \
+ > $@ || (rm -f $@; false)
+ @if test -s $@ ; then : else rm $@ ; fi
+
+include Makefile.dep
Added: cs/cigma/trunk/examples/laplace1/README
===================================================================
--- cs/cigma/trunk/examples/laplace1/README (rev 0)
+++ cs/cigma/trunk/examples/laplace1/README 2009-01-13 11:11:59 UTC (rev 13857)
@@ -0,0 +1 @@
+To run this example, you need to install the deal.II library
Added: cs/cigma/trunk/examples/laplace1/laplace.cc
===================================================================
--- cs/cigma/trunk/examples/laplace1/laplace.cc (rev 0)
+++ cs/cigma/trunk/examples/laplace1/laplace.cc 2009-01-13 11:11:59 UTC (rev 13857)
@@ -0,0 +1,267 @@
+//
+// Modified step-4.cc example from deal.II source
+//
+#include <grid/tria.h>
+#include <dofs/dof_handler.h>
+#include <grid/grid_generator.h>
+#include <grid/tria_accessor.h>
+#include <grid/tria_iterator.h>
+#include <dofs/dof_accessor.h>
+#include <fe/fe_q.h>
+#include <dofs/dof_tools.h>
+#include <fe/fe_values.h>
+#include <base/quadrature_lib.h>
+#include <base/function.h>
+#include <numerics/vectors.h>
+#include <numerics/matrices.h>
+#include <lac/vector.h>
+#include <lac/full_matrix.h>
+#include <lac/sparse_matrix.h>
+#include <lac/solver_cg.h>
+#include <lac/precondition.h>
+
+#include <numerics/data_out.h>
+#include <iostream>
+#include <fstream>
+#include <sstream>
+
+#include <base/logstream.h>
+
+using namespace std;
+using namespace dealii;
+
+// ----------------------------------------------------------------------------
+// Function representing the right hand side term in our equation.
+//
+template <int dim>
+class RightHandSide : public Function<dim>
+{
+public:
+ RightHandSide() : Function<dim>() {}
+ virtual double value(const Point<dim> &p, const unsigned int component=0) const;
+};
+
+template <int dim>
+double RightHandSide<dim>::value(const Point<dim> &p, const unsigned int /*component*/) const
+{
+ double return_value = 0;
+ for (unsigned int i = 0; i < dim; i++)
+ {
+ return_value += 4 * std::pow(p(i), 4);
+ }
+ return return_value;
+}
+
+
+// ----------------------------------------------------------------------------
+// Function for calculating function at boundary.
+//
+template <int dim>
+class BoundaryValues : public Function<dim>
+{
+public:
+ BoundaryValues() : Function<dim>() {}
+ virtual double value(const Point<dim> &p, const unsigned int component=0) const;
+};
+
+template <int dim>
+double BoundaryValues<dim>::value(const Point<dim> &p, const unsigned int /*component*/) const
+{
+ return p.square();
+}
+
+
+// ----------------------------------------------------------------------------
+// This class solves the Laplace problem
+//
+template <int dim>
+class LaplaceProblem
+{
+public:
+ LaplaceProblem();
+ void run();
+
+private:
+ void make_grid_and_dofs();
+ void assemble_system();
+ void solve();
+ void output_results() const;
+
+private:
+ Triangulation<dim> triangulation;
+ FE_Q<dim> fe;
+ DoFHandler<dim> dof_handler;
+
+ SparsityPattern sparsity_pattern;
+ SparseMatrix<double> system_matrix;
+
+ Vector<double> solution;
+ Vector<double> system_rhs;
+};
+
+template <int dim>
+LaplaceProblem<dim>::LaplaceProblem() : fe(1), dof_handler(triangulation)
+{
+}
+
+template <int dim>
+void LaplaceProblem<dim>::make_grid_and_dofs()
+{
+ GridGenerator::hyper_cube(triangulation, -1, +1);
+ triangulation.refine_global(3);
+
+ cout << "\tNumber of active cells: " << triangulation.n_active_cells() << endl;
+ cout << "\tTotal number of cells: " << triangulation.n_cells() << endl;
+
+ dof_handler.distribute_dofs(fe);
+
+ cout << "\tNumber of degrees of freedom: " << dof_handler.n_dofs() << endl;
+
+ sparsity_pattern.reinit(
+ dof_handler.n_dofs(),
+ dof_handler.n_dofs(),
+ dof_handler.max_couplings_between_dofs());
+
+ DoFTools::make_sparsity_pattern(dof_handler, sparsity_pattern);
+ sparsity_pattern.compress();
+
+ system_matrix.reinit(sparsity_pattern);
+
+ solution.reinit(dof_handler.n_dofs());
+ system_rhs.reinit(dof_handler.n_dofs());
+}
+
+template <int dim>
+void LaplaceProblem<dim>::assemble_system()
+{
+ QGauss<dim> quadrature_formula(2);
+ const RightHandSide<dim> right_hand_side;
+ FEValues<dim> fe_values(fe,
+ quadrature_formula,
+ update_values |
+ update_gradients |
+ update_quadrature_points |
+ update_JxW_values);
+
+ const unsigned int dofs_per_cell = fe.dofs_per_cell;
+ const unsigned int n_q_points = quadrature_formula.size();
+
+ FullMatrix<double> cell_matrix(dofs_per_cell, dofs_per_cell);
+ Vector<double> cell_rhs(dofs_per_cell);
+
+ std::vector<unsigned int> local_dof_indices(dofs_per_cell);
+
+ typename DoFHandler<dim>::active_cell_iterator
+ cell = dof_handler.begin_active(),
+ endc = dof_handler.end();
+
+ for (; cell != endc; ++cell)
+ {
+ fe_values.reinit(cell);
+ cell_matrix = 0;
+ cell_rhs = 0;
+
+ for (unsigned int q_point = 0; q_point < n_q_points; q_point++)
+ {
+ for (unsigned int i = 0; i < dofs_per_cell; i++)
+ {
+ for (unsigned int j = 0; j < dofs_per_cell; j++)
+ {
+ cell_matrix(i,j) +=
+ fe_values.shape_grad(i, q_point) *
+ fe_values.shape_grad(j, q_point) *
+ fe_values.JxW(q_point);
+ }
+ cell_rhs(i) +=
+ fe_values.shape_value(i, q_point) *
+ right_hand_side.value(fe_values.quadrature_point(q_point)) *
+ fe_values.JxW(q_point);
+ }
+
+ cell->get_dof_indices(local_dof_indices);
+
+ for(unsigned int i = 0; i < dofs_per_cell; i++)
+ {
+ for (unsigned int j = 0; j < dofs_per_cell; j++)
+ {
+ system_matrix.add(local_dof_indices[i], local_dof_indices[j], cell_matrix(i,j));
+ }
+ system_rhs(local_dof_indices[i]) += cell_rhs(i);
+ }
+ }
+ }
+
+ std::map<unsigned int, double> boundary_values;
+
+ VectorTools::interpolate_boundary_values(
+ dof_handler,
+ 0,
+ BoundaryValues<dim>(),
+ boundary_values);
+
+ MatrixTools::apply_boundary_values(
+ boundary_values,
+ system_matrix,
+ solution,
+ system_rhs);
+}
+
+
+template <int dim>
+void LaplaceProblem<dim>::solve()
+{
+ SolverControl solver_control(1000, 1e-12);
+ SolverCG<> cg(solver_control);
+ cg.solve(system_matrix, solution, system_rhs, PreconditionIdentity());
+
+ cout << "\tNeeded " << solver_control.last_step()
+ << " CG iterations to obtain convergence."
+ << endl;
+}
+
+template <int dim>
+void LaplaceProblem<dim>::output_results() const
+{
+ DataOut<dim> data_out;
+ data_out.attach_dof_handler(dof_handler);
+ data_out.add_data_vector(solution, "solution");
+ data_out.build_patches();
+
+ std::ostringstream name;
+ name << "solution-" << dim << "d" << ".vtk" << std::ends;
+ std::ofstream output(name.str().c_str());
+ data_out.write_vtk(output);
+}
+
+template <int dim>
+void LaplaceProblem<dim>::run()
+{
+ cout << "Solving problem in " << dim << " space dimensions." << endl;
+
+ make_grid_and_dofs();
+ assemble_system();
+ solve();
+ output_results();
+}
+
+// ----------------------------------------------------------------------------
+// The main function
+
+int main(void)
+{
+ deallog.depth_console(0);
+
+ if (true)
+ {
+ LaplaceProblem<2> laplace_problem_in_2d;
+ laplace_problem_in_2d.run();
+ }
+
+ if (true)
+ {
+ LaplaceProblem<3> laplace_problem_in_3d;
+ laplace_problem_in_3d.run();
+ }
+
+ return 0;
+}
Added: cs/cigma/trunk/examples/laplace2/Makefile
===================================================================
--- cs/cigma/trunk/examples/laplace2/Makefile (rev 0)
+++ cs/cigma/trunk/examples/laplace2/Makefile 2009-01-13 11:11:59 UTC (rev 13857)
@@ -0,0 +1,61 @@
+###############################################################################
+# -*- Makefile -*-
+
+# The location of the mesh library
+meshdir := $(HOME)/dev/libmesh
+
+# Include the library options determined by configure.
+# This will set the variables INCLUDE and LIBS that we
+# will need to build and link with the library
+include $(meshdir)/Make.common
+
+target := ./laplace-$(METHOD)
+
+# source files
+srcfiles := $(wildcard *.cpp)
+
+# object files
+objects := $(patsubst %.cpp, %.$(obj-suffix), $(srcfiles))
+
+
+all:: $(target)
+
+# Production rules
+$(target): $(objects)
+ @echo Linking "$@"...
+ $(libmesh_CXX) $(libmesh_CXXFLAGS) $(objects) -o $@ $(libmesh_LIBS) $(libmesh_LDFLAGS)
+
+# Rules for running examples
+run: run1 run2 run3
+
+run1: $(target)
+ $(LIBMESHRUN) $(target) -d 1 -n 20 $(LIBMESHOPTIONS)
+
+run2: $(target)
+ $(LIBMESHRUN) $(target) -d 1 -n 15 $(LIBMESHOPTIONS)
+
+run3: $(target)
+ $(LIBMESHRUN) $(target) -d 1 -n 6 $(LIBMESHOPTIONS)
+
+# Useful rules
+clean:
+ rm -f $(objects) *~
+
+clobber:
+ @$(MAKE) clean
+ rm -f $(target) out_*
+
+distclean:
+ @$(MAKE) clobber
+ rm -f *.o *.g.o *.pg.o
+
+.PHONY: clean clobber distclean
+
+
+# include the dependency list
+include Makefile.dep
+
+Makefile.dep:
+ @$(perl) $(meshdir)/contrib/bin/make_dependencies.pl -I. $(foreach i, $(wildcard $(meshdir)/include/*), -I$(i)) "-S\$$(obj-suffix)" $(srcfiles) > $@
+ @echo Updated .depend
+
Added: cs/cigma/trunk/examples/laplace2/README
===================================================================
--- cs/cigma/trunk/examples/laplace2/README (rev 0)
+++ cs/cigma/trunk/examples/laplace2/README 2009-01-13 11:11:59 UTC (rev 13857)
@@ -0,0 +1,5 @@
+To run this example, you need to install libMesh. The libMesh library is
+a C++ framework for the numerical simulation of partial differential equations.
+
+On a Debian system, it is sufficient to install the libmesh-dev package.
+
Added: cs/cigma/trunk/examples/laplace2/exact_solution.C
===================================================================
--- cs/cigma/trunk/examples/laplace2/exact_solution.C (rev 0)
+++ cs/cigma/trunk/examples/laplace2/exact_solution.C 2009-01-13 11:11:59 UTC (rev 13857)
@@ -0,0 +1,24 @@
+#include "exact_solution.h"
+#include <cmath>
+
+/**
+ * This is the exact solution we are trying to obtain. We solve
+ *
+ * - (u_xx + u_yy) = f
+ *
+ * and take a finite difference approximation using this function
+ * to get f. This is the well-known "method of manufactured solutions".
+ */
+
+Real exact_solution(const Real x, const Real y, const Real z=0.0)
+{
+ static const Real pi = acos(-1.0);
+ //return cos(0.5*pi*x)*sin(0.5*pi*y)*cos(0.5*pi*z);
+ return 4*(x*x*x*x + y*y*y*y + z*z*z*z);
+}
+
+Real exact_boundary(const Real x, const Real y=0.0, const Real z=0.0)
+{
+ return x*x + y*y + z*z;
+}
+
Added: cs/cigma/trunk/examples/laplace2/exact_solution.H
===================================================================
--- cs/cigma/trunk/examples/laplace2/exact_solution.H (rev 0)
+++ cs/cigma/trunk/examples/laplace2/exact_solution.H 2009-01-13 11:11:59 UTC (rev 13857)
@@ -0,0 +1,9 @@
+#ifndef EXACT_SOLUTION_H
+#define EXACT_SOLUTION_H
+
+#include "libmesh_common.h"
+
+Real exact_solution(const Real x, const Real y, const Real z=0.0);
+Real exact_boundary(const Real x, const Real y=0.0, const Real z=0.0);
+
+#endif
Added: cs/cigma/trunk/examples/laplace2/laplace.C
===================================================================
--- cs/cigma/trunk/examples/laplace2/laplace.C (rev 0)
+++ cs/cigma/trunk/examples/laplace2/laplace.C 2009-01-13 11:11:59 UTC (rev 13857)
@@ -0,0 +1,369 @@
+// Modified version of ex4.C example in libMesh source code
+//
+
+#include <iostream>
+#include <algorithm>
+#include <sstream>
+#include <cmath>
+
+#include "libmesh.h"
+#include "mesh.h"
+#include "mesh_generation.h"
+#include "vtk_io.h"
+#include "ucd_io.h"
+#include "gmv_io.h"
+#include "gnuplot_io.h"
+#include "linear_implicit_system.h"
+#include "equation_systems.h"
+
+#include "fe.h"
+#include "quadrature_gauss.h"
+#include "dof_map.h"
+
+// Useful datatypes for finite element matrix and vector components
+#include "sparse_matrix.h"
+#include "numeric_vector.h"
+#include "dense_matrix.h"
+#include "dense_vector.h"
+
+// defines PerfLog
+#include "perf_log.h"
+
+// definition of a geometric element
+#include "elem.h"
+
+// misc
+#include "string_to_enum.h"
+#include "getpot.h"
+
+// Exact solution
+#include "exact_solution.h"
+
+// Assemble the linear system for Poisson problem
+void assemble_poisson(EquationSystems& es, const std::string& system_name)
+{
+ libmesh_assert(system_name == "Poisson");
+
+ PerfLog perf_log("Matrix Assembly");
+
+ const MeshBase& mesh = es.get_mesh();
+ const unsigned int dim = mesh.mesh_dimension();
+
+ LinearImplicitSystem& system = es.get_system<LinearImplicitSystem>("Poisson");
+ const DofMap& dof_map = system.get_dof_map();
+
+ FEType fe_type = dof_map.variable_type(0);
+
+ AutoPtr<FEBase> fe(FEBase::build(dim, fe_type));
+ QGauss qrule(dim, FIFTH);
+ fe->attach_quadrature_rule(&qrule);
+
+ AutoPtr<FEBase> fe_face(FEBase::build(dim, fe_type));
+ QGauss qface(dim-1, FIFTH);
+ fe_face->attach_quadrature_rule(&qface);
+
+ // Some typedefs
+ typedef std::vector<Real> VectorOfReal;
+ typedef std::vector<RealGradient> VectorOfRealGradient;
+ typedef std::vector<Point> VectorOfPoint;
+
+ // Element jacobian * quadrature weight at each integration point
+ const VectorOfReal& JxW = fe->get_JxW();
+
+ // Physical XYZ locations of the quadrature points on the element
+ const VectorOfPoint& q_point = fe->get_xyz();
+
+ // Element shape functions evaluated at the quadrature points
+ const std::vector<VectorOfReal>& phi = fe->get_phi();
+
+ // Element shape function gradients evalauted at the quadrature points;
+ const std::vector<VectorOfRealGradient>& dphi = fe->get_dphi();
+
+ // Data structures for element matrix and right-hand-side vector contribution.
+ DenseMatrix<Number> Ke;
+ DenseVector<Number> Fe;
+
+ // Hold the degree of freedom indices for the element
+ std::vector<unsigned int> dof_indices;
+
+ // Loop over all elements in the mesh
+ MeshBase::const_element_iterator el = mesh.local_elements_begin();
+ const MeshBase::const_element_iterator end_el = mesh.local_elements_end();
+
+ for ( ; el != end_el; el++)
+ {
+ // Start logging the shape function initialization
+ perf_log.push("elem init");
+ const Elem* elem = *el;
+
+ // Get the degree of freedom indices for the current element
+ dof_map.dof_indices(elem, dof_indices);
+
+ // Compute the element-specific data for the current element
+ fe->reinit(elem);
+ Ke.resize(dof_indices.size(), dof_indices.size());
+ Fe.resize(dof_indices.size());
+ perf_log.pop("elem init");
+
+
+ // Start logging the element matrix computation
+ perf_log.push("Ke");
+ for (unsigned int qp = 0; qp < qrule.n_points(); qp++)
+ {
+ for (unsigned int i = 0; i < phi.size(); i++)
+ {
+ for (unsigned int j = 0; j < phi.size(); j++)
+ {
+ Ke(i,j) += JxW[qp] * dphi[i][qp] * dphi[j][qp];
+ }
+ }
+ }
+ perf_log.pop("Ke");
+
+ // Start logging the right-hand-side computation
+ perf_log.push("Fe");
+ for (unsigned int qp = 0; qp < qrule.n_points(); qp++)
+ {
+ const Real x = q_point[qp](0);
+ const Real y = q_point[qp](1);
+ const Real z = q_point[qp](2);
+ const Real eps = 1.0e-3;
+
+ const Real uxx =
+ (exact_solution(x-eps, y, z) +
+ exact_solution(x+eps, y, z) +
+ -2.0*exact_solution(x,y,z))/eps/eps;
+
+ const Real uyy =
+ (exact_solution(x, y-eps, z) +
+ exact_solution(x, y+eps, z) +
+ -2.0*exact_solution(x,y,z))/eps/eps;
+
+ const Real uzz =
+ (exact_solution(x, y, z-eps) +
+ exact_solution(x, y, z+eps) +
+ -2.0*exact_solution(x,y,z))/eps/eps;
+
+ Real fxy;
+ if (dim == 1)
+ {
+ // In 1d, compute the rhs by differentiating
+ // the exact solution twice
+ const Real pi = libMesh::pi;
+ fxy = 0.25*pi*pi*sin(0.5*pi*x);
+ }
+ else if (dim == 2)
+ {
+ fxy = -(uxx + uyy);
+ }
+ else
+ {
+ fxy = -(uxx + uyy + uzz);
+ }
+
+ }
+ perf_log.pop("Fe");
+
+ // Interior element integration is completed.
+ {
+ // Start logging the boundary condition computation
+ perf_log.push("BCs");
+ for (unsigned int side = 0; side < elem->n_sides(); side++)
+ {
+ if (elem->neighbor(side) == NULL)
+ {
+ const Real penalty = 1e10;
+ const std::vector<VectorOfReal>& phi_face = fe_face->get_phi();
+ const VectorOfReal& JxW_face = fe_face->get_JxW();
+ const VectorOfPoint& qface_point = fe_face->get_xyz();
+
+ fe_face->reinit(elem, side);
+
+ for (unsigned int qp = 0; qp < qface.n_points(); qp++)
+ {
+ // The location on the boundary of the current face quadrature point
+ const Real xf = qface_point[qp](0);
+ const Real yf = qface_point[qp](1);
+ const Real zf = qface_point[qp](2);
+
+ // The boundary value
+ const Real value = exact_boundary(xf, yf, zf);
+
+ // Matrix contribution of the L2 projection
+ for (unsigned int i = 0; i < phi_face.size(); i++)
+ {
+ for (unsigned int j = 0; j < phi_face.size(); j++)
+ {
+ Ke(i,j) += JxW_face[qp] * penalty * phi_face[i][qp] * phi_face[j][qp];
+ }
+ }
+
+ // The right-hand-side contribution of the L2 projection
+ for (unsigned int i = 0; i < phi_face.size(); i++)
+ {
+ Fe(i) += JxW_face[qp] * penalty * value * phi_face[i][qp];
+ }
+ }
+ }
+ }
+ perf_log.pop("BCs");
+ }
+
+ // The element matrix and right-hand-side are now built for this element
+ perf_log.push("matrix insertion");
+
+ system.matrix->add_matrix(Ke, dof_indices);
+ system.rhs->add_vector(Fe, dof_indices);
+
+ perf_log.pop("matrix insertion");
+ }
+
+ return;
+}
+
+
+using namespace std;
+
+int main(int argc, char *argv[])
+{
+ LibMeshInit init(argc, argv);
+ PerfLog perf_main("Main Program");
+
+ GetPot command_line(argc, argv);
+
+ if (argc < 3)
+ {
+ if (libMesh::processor_id() == 0)
+ {
+ cerr << "Usage:\n\t" << argv[0] << " -d [ 2 | 3] -n 15" << endl;
+ libmesh_error();
+ }
+ }
+ else
+ {
+ cout << "Running " << argv[0];
+ for (int i = 1; i < argc; i++)
+ cout << " " argv[i];
+ cout << endl << endl;
+ }
+
+ // Read output filename from command line
+ string outfile = "solution.pvtu";
+ if (command_line.search(1, "-t"))
+ {
+ outfile = command_line.next(outfile);
+ }
+
+ // Read problem dimension from command line
+ int dim = 2;
+ if (command_line.search(1, "-d"))
+ {
+ dim = command_line.next(dim);
+ }
+
+ // Read number of elements from command line
+ int ps = 15;
+ if (command_line.search(1, "-n"))
+ {
+ ps = command_line.next(ps);
+ }
+
+ // Read FE order from command line
+ string order = "FIRST";
+ //string order = "SECOND";
+ if (command_line.search(2, "-Order", "-o"))
+ {
+ order = command_line.next(order);
+ }
+
+ // Read FE family from command line
+ string family = "LAGRANGE";
+ if (command_line.search(2, "-FEFamily", "-f"))
+ {
+ family = command_line.next(family);
+ }
+
+ // Don't use a discontinuous basis
+ if ((family == "MONOMIAL") || (family = "XYZ"))
+ {
+ cout << "This example currently requires a C^0 (or higher) FE basis." << endl;
+ lib_mesh_error();
+ }
+
+ // Create a mesh with user-defined dimension
+ Mesh mesh(dim);
+
+ if ((family == "LAGRANGE") && (order == "FIRST"))
+ {
+ MeshTools::Generation::build_cube(
+ mesh, ps, ps, ps,
+ -1.0, 1.0,
+ -1.0, 1.0,
+ -1.0, 1.0,
+ (dim == 1) ? EDGE2 : ((dim == 2) ? QUAD4 : HEX8)
+ );
+ }
+ else
+ {
+ MeshTools::Generation::build_cube(
+ mesh,
+ ps, ps, ps,
+ -1.0, 1.0,
+ -1.0, 1.0,
+ -1.0, 1.0,
+ (dim == 1) ? EDGE3 : ((dim == 2) ? QUAD9 : HEX27)
+ );
+ }
+
+ mesh.print_info();
+
+ // Create an equation systems object
+ EquationSystems equation_systems(mesh);
+
+ LinearImplicitSystem& system = equation_systems.add_system<LinearImplicitSystem>("Poisson");
+
+ system.add_variable(
+ "u",
+ Utility::string_to_enum<Order>(order),
+ Utility::string_to_enum<FEFamily>(family));
+
+ system.attach_assemble_function(assemble_poisson);
+
+ equation_systems.init();
+
+ equation_systems.print_info();
+ mesh.print_info();
+
+ equation_systems.get_system("Poisson").solve();
+
+ if (dim == 1)
+ {
+ GnuPlotIO plot(mesh, "Laplace Example, 1D", GnuPlotIO::GRID_ON);
+ plot.write_equation_systems("solution1", equation_systems);
+ }
+ else
+ {
+ string ext = ".vtk";
+
+ ostringstream stream;
+ stream << "solution-" << dim << ext << std::ends;
+
+ const char *filename = stream.str().c_str();
+
+ if (ext == ".vtk")
+ {
+ VTKIO(mesh).write_equation_systems(filename, equation_systems);
+ }
+ else if (ext == ".gmv")
+ {
+ GMVIO(mesh).write_equation_systems(filename, equation_systems);
+ }
+ else if (ext == ".ucd")
+ {
+ UCDIO(mesh).write_equation_systems(filename, equation_systems);
+ }
+
+ }
+
+ return 0;
+}
More information about the CIG-COMMITS
mailing list