On Mon, 29 Oct 2012, Kirk, Benjamin (JSC-EG311) wrote:
I don't think there are any examples that currently show this - I'll
try to add one after the next release.
See the attached ~300 lines of code. It reads a mesh and solution, L2
projects and then plots all gradient components of all solution
variables. It needs to be changed to have more than just
Tecplot-or-GMV output, and as you may have guessed this implies I
haven't tested it in a couple years, but other than that it ought to
be ready to drop in as fem_system_ex2.
---
Roy
// Open the mesh and solution file named in command line arguments,
// L2 project their gradient components,
// and plot them out to the given filename.
#include "libmesh.h"
#include "dof_map.h"
#include "equation_systems.h"
#include "fe_base.h"
#include "fem_context.h"
#include "fem_system.h"
#include "mesh.h"
#include "quadrature.h"
#include "system.h"
#include "gmv_io.h"
#include "tecplot_io.h"
class GradientSystem : public FEMSystem
{
public:
// Constructor
GradientSystem(EquationSystems& es,
const std::string& name,
const unsigned int number)
: FEMSystem(es, name, number) {}
System *original_system;
protected:
// System initialization
virtual void init_data ();
// Context initialization
virtual void init_context (DiffContext &context);
// Element residual and jacobian calculations
virtual bool element_constraint (bool request_jacobian,
DiffContext &context);
};
void GradientSystem::init_data ()
{
unsigned int n_original_vars = original_system->n_vars();
for (unsigned int v=0; v != n_original_vars; ++v)
{
std::string varname = original_system->variable_name(v);
this->add_variable(varname+"_x",
original_system->variable_type(v).order,
original_system->variable_type(v).family);
#if LIBMESH_DIM > 1
this->add_variable(varname+"_y",
original_system->variable_type(v).order,
original_system->variable_type(v).family);
#endif
#if LIBMESH_DIM > 2
this->add_variable(varname+"_z",
original_system->variable_type(v).order,
original_system->variable_type(v).family);
#endif
}
// Do the parent's initialization after variables are defined
FEMSystem::init_data();
}
void GradientSystem::init_context(DiffContext &context)
{
FEMContext &c = libmesh_cast_ref<FEMContext&>(context);
unsigned int n_original_vars = original_system->n_vars();
// Make sure we have requested all the data we need to build each
// linear system. Since we're reusing FEs of the same FEType, we
// don't need to do separate _x,_y,_z requests.
for (unsigned int v=0; v != n_original_vars; ++v)
{
c.element_fe_var[v]->get_JxW();
c.element_fe_var[v]->get_phi();
c.element_fe_var[v]->get_dphi();
}
FEMSystem::init_context(context);
}
bool GradientSystem::element_constraint (bool request_jacobian,
DiffContext &context)
{
FEMContext &c = libmesh_cast_ref<FEMContext&>(context);
unsigned int n_original_vars = original_system->n_vars();
// We assemble a separate linear system for each gradient direction
// of each variable. Solving this fully-coupled is unnecessary but
// doesn't hurt much.
for (unsigned int v = 0; v != n_original_vars; ++v)
{
// We have LIBMESH_DIM gradient variables for each original
// variable
unsigned int vx = v * LIBMESH_DIM;
#if LIBMESH_DIM > 1
unsigned int vy = v * LIBMESH_DIM+1;
#endif
#if LIBMESH_DIM > 2
unsigned int vz = v * LIBMESH_DIM+2;
#endif
// Element Jacobian * quadrature weights for interior integration
const std::vector<Real> &JxW = c.element_fe_var[v]->get_JxW();
// For test function values
const std::vector<std::vector<Real> > &phi =
c.element_fe_var[v]->get_phi();
// For target gradient values
const std::vector<std::vector<RealGradient> > &dphi =
c.element_fe_var[v]->get_dphi();
// The target dof indices
std::vector<unsigned int> original_dof_indices;
original_system->get_dof_map().dof_indices(c.elem, original_dof_indices,
v);
// The number of local degrees of freedom
const unsigned int n_u_dofs = c.dof_indices_var[vx].size();
// We're using the same basis for the original system and for
// each gradient component
libmesh_assert(n_u_dofs == original_dof_indices.size());
#if LIBMESH_DIM > 1
libmesh_assert(n_u_dofs == c.dof_indices_var[vy].size());
#endif
#if LIBMESH_DIM > 2
libmesh_assert(n_u_dofs == c.dof_indices_var[vz].size());
#endif
// The subvectors and submatrices we need to fill:
DenseSubMatrix<Number> &Kxx = *c.elem_subjacobians[vx][vx];
DenseSubVector<Number> &Fx = *c.elem_subresiduals[vx];
#if LIBMESH_DIM > 1
DenseSubMatrix<Number> &Kyy = *c.elem_subjacobians[vy][vy];
DenseSubVector<Number> &Fy = *c.elem_subresiduals[vy];
#endif
#if LIBMESH_DIM > 2
DenseSubMatrix<Number> &Kzz = *c.elem_subjacobians[vz][vz];
DenseSubVector<Number> &Fz = *c.elem_subresiduals[vz];
#endif
// Now we will build the element Jacobian and residual.
// Constructing the residual requires the solution and its
// gradient from the previous timestep. This must be
// calculated at each quadrature point by summing the
// solution degree-of-freedom values by the appropriate
// weight functions.
const unsigned int n_qpoints = c.element_qrule->n_points();
for (unsigned int qp=0; qp != n_qpoints; qp++)
{
// Compute the gradient we want to project
Gradient g = 0.;
for (unsigned int l=0; l != n_u_dofs; l++)
{
g.add_scaled(dphi[l][qp],
original_system->current_solution(original_dof_indices[l]));
}
// Compute the solution at the Newton iterate; we don't want
// the solution to drift if the user uses too weak a linear
// tolerance
const Number ux_iter = c.interior_value(vx, qp);
#if LIBMESH_DIM > 1
const Number uy_iter = c.interior_value(vy, qp);
#endif
#if LIBMESH_DIM > 2
const Number uz_iter = c.interior_value(vz, qp);
#endif
for (unsigned int i=0; i != n_u_dofs; i++)
{
Fx(i) += JxW[qp] * (ux_iter - g(0)) * phi[i][qp];
#if LIBMESH_DIM > 1
Fy(i) += JxW[qp] * (uy_iter - g(1)) * phi[i][qp];
#endif
#if LIBMESH_DIM > 2
Fz(i) += JxW[qp] * (uz_iter - g(2)) * phi[i][qp];
#endif
}
if (request_jacobian)
for (unsigned int i=0; i != n_u_dofs; i++)
for (unsigned int j=0; j != n_u_dofs; ++j)
{
Kxx(i,j) += JxW[qp] * (phi[i][qp] * phi[j][qp]);
#if LIBMESH_DIM > 1
Kyy(i,j) += JxW[qp] * (phi[i][qp] * phi[j][qp]);
#endif
#if LIBMESH_DIM > 2
Kzz(i,j) += JxW[qp] * (phi[i][qp] * phi[j][qp]);
#endif
}
} // end of the quadrature point qp-loop
} // end of the variable v-loop
return request_jacobian;
}
int main(int argc, char** argv)
{
LibMeshInit init(argc, argv);
Mesh mesh;
EquationSystems es(mesh);
std::cout << "Usage: " << argv[0]
<< " inputmesh inputsolution outputsolution" << std::endl;
START_LOG("mesh.read()", "main");
mesh.read(argv[1]);
std::cout << "Loaded mesh " << argv[1] << std::endl;
STOP_LOG("mesh.read()", "main");
START_LOG("es.read()", "main");
std::string solnname = argv[2];
libmesh_assert(solnname.find(".xdr") != std::string::npos ||
solnname.find(".xda") != std::string::npos);
const bool binary = (solnname.find(".xdr") != std::string::npos);
es.read(argv[2], binary ? libMeshEnums::DECODE : libMeshEnums::READ,
EquationSystems::READ_HEADER |
EquationSystems::READ_DATA |
EquationSystems::READ_ADDITIONAL_DATA);
std::cout << "Loaded " << (binary ? "binary" : "ascii") << " solution " <<
argv[2] << std::endl;
STOP_LOG("es.read()", "main");
const unsigned int n_systems = es.n_systems();
// We'll use a copy of the mesh and a new EquationSystems so as to
// leave the original Systems out of our writes
AutoPtr<MeshBase> newmesh = mesh.clone();
EquationSystems newes(*newmesh);
for (unsigned int s = 0; s != n_systems; ++s)
{
System &orig = es.get_system(s);
std::string sysname = orig.name();
GradientSystem &gradsys = newes.add_system<GradientSystem>(sysname);
gradsys.original_system = &orig;
}
newes.solve();
START_LOG("write_equation_systems()", "main");
std::string outputname(argv[3]);
if (outputname.find(".gmv") != std::string::npos)
GMVIO(*newmesh).write_equation_systems (outputname, newes);
else if (outputname.find(".plt") != std::string::npos)
TecplotIO(*newmesh).write_equation_systems (outputname, es);
STOP_LOG("write_equation_systems()", "main");
std::cout << "Wrote output " << argv[3] << std::endl;
}
------------------------------------------------------------------------------
The Windows 8 Center - In partnership with Sourceforge
Your idea - your app - 30 days.
Get started!
http://windows8center.sourceforge.net/
what-html-developers-need-to-know-about-coding-windows-8-metro-style-apps/
_______________________________________________
Libmesh-devel mailing list
Libmesh-devel@lists.sourceforge.net
https://lists.sourceforge.net/lists/listinfo/libmesh-devel
