libMesh::ProjectSolution Class Reference
Public Member Functions | |
| ProjectSolution (const System &system_in, FunctionBase< Number > *f_in, FunctionBase< Gradient > *g_in, const Parameters ¶meters_in, NumericVector< Number > &new_v_in) | |
| ProjectSolution (const ProjectSolution &in) | |
| void | operator() (const ConstElemRange &range) const |
Private Attributes | |
| const System & | system |
| UniquePtr< FunctionBase< Number > > | f |
| UniquePtr< FunctionBase < Gradient > > | g |
| const Parameters & | parameters |
| NumericVector< Number > & | new_vector |
Detailed Description
This class implements projecting an arbitrary function to the current mesh. This may be exectued in parallel on multiple threads.
Definition at line 110 of file system_projection.C.
Constructor & Destructor Documentation
| libMesh::ProjectSolution::ProjectSolution | ( | const System & | system_in, |
| FunctionBase< Number > * | f_in, | ||
| FunctionBase< Gradient > * | g_in, | ||
| const Parameters & | parameters_in, | ||
| NumericVector< Number > & | new_v_in | ||
| ) | [inline] |
Definition at line 121 of file system_projection.C.
References f, g, and libMesh::libmesh_assert().
:
system(system_in),
f(f_in ? f_in->clone() : UniquePtr<FunctionBase<Number> >()),
g(g_in ? g_in->clone() : UniquePtr<FunctionBase<Gradient> >()),
parameters(parameters_in),
new_vector(new_v_in)
{
libmesh_assert(f.get());
f->init();
if (g.get())
g->init();
}
| libMesh::ProjectSolution::ProjectSolution | ( | const ProjectSolution & | in | ) | [inline] |
Definition at line 138 of file system_projection.C.
References f, g, and libMesh::libmesh_assert().
:
system(in.system),
f(in.f.get() ? in.f->clone() : UniquePtr<FunctionBase<Number> >()),
g(in.g.get() ? in.g->clone() : UniquePtr<FunctionBase<Gradient> >()),
parameters(in.parameters),
new_vector(in.new_vector)
{
libmesh_assert(f.get());
f->init();
if (g.get())
g->init();
}
Member Function Documentation
| void libMesh::ProjectSolution::operator() | ( | const ConstElemRange & | range | ) | const |
This method projects an arbitrary solution to the current mesh. The input function f gives the arbitrary solution, while the new_vector (which should already be correctly sized) gives the solution (to be computed) on the current mesh.
Definition at line 1270 of file system_projection.C.
References std::abs(), libMesh::Variable::active_on_subdomain(), libMesh::StoredRange< iterator_type, object_type >::begin(), libMesh::FEGenericBase< OutputType >::build(), libMesh::C_ONE, libMesh::C_ZERO, libMesh::DenseMatrix< T >::cholesky_solve(), libMesh::FunctionBase< Output >::component(), libMesh::FEType::default_quadrature_rule(), libMesh::DISCONTINUOUS, libMesh::DofMap::dof_indices(), libMesh::FEInterface::dofs_on_edge(), libMesh::FEInterface::dofs_on_side(), libMesh::StoredRange< iterator_type, object_type >::end(), libMesh::FEType::family, libMesh::NumericVector< T >::first_local_index(), libMesh::System::get_dof_map(), libMesh::System::get_mesh(), libMesh::HERMITE, libMesh::Elem::is_vertex(), libMesh::NumericVector< T >::last_local_index(), libMesh::libmesh_assert(), libMesh::MeshBase::mesh_dimension(), libMesh::FEInterface::n_dofs_at_node(), libMesh::Elem::n_edges(), n_nodes, libMesh::Elem::n_nodes(), libMesh::Elem::n_sides(), libMesh::System::n_vars(), libMesh::Elem::point(), libMesh::DenseVector< T >::resize(), libMesh::DenseMatrix< T >::resize(), libMesh::SCALAR, libMesh::NumericVector< T >::set(), libMesh::Threads::spin_mtx, libMesh::Elem::subdomain_id(), system, libMesh::System::time, libMesh::TOLERANCE, libMesh::Variable::type(), libMesh::Elem::type(), libMesh::DofMap::variable(), libMesh::System::variable_scalar_number(), libMesh::DenseMatrix< T >::zero(), and libMesh::DenseVector< T >::zero().
{
// We need data to project
libmesh_assert(f.get());
// The number of variables in this system
const unsigned int n_variables = system.n_vars();
// The dimensionality of the current mesh
const unsigned int dim = system.get_mesh().mesh_dimension();
// The DofMap for this system
const DofMap& dof_map = system.get_dof_map();
// The element matrix and RHS for projections.
// Note that Ke is always real-valued, whereas
// Fe may be complex valued if complex number
// support is enabled
DenseMatrix<Real> Ke;
DenseVector<Number> Fe;
// The new element coefficients
DenseVector<Number> Ue;
// Loop over all the variables in the system
for (unsigned int var=0; var<n_variables; var++)
{
const Variable& variable = dof_map.variable(var);
const FEType& fe_type = variable.type();
if (fe_type.family == SCALAR)
continue;
const unsigned int var_component =
system.variable_scalar_number(var, 0);
// Get FE objects of the appropriate type
UniquePtr<FEBase> fe (FEBase::build(dim, fe_type));
// Prepare variables for projection
UniquePtr<QBase> qrule (fe_type.default_quadrature_rule(dim));
UniquePtr<QBase> qedgerule (fe_type.default_quadrature_rule(1));
UniquePtr<QBase> qsiderule (fe_type.default_quadrature_rule(dim-1));
// The values of the shape functions at the quadrature
// points
const std::vector<std::vector<Real> >& phi = fe->get_phi();
// The gradients of the shape functions at the quadrature
// points on the child element.
const std::vector<std::vector<RealGradient> > *dphi = NULL;
const FEContinuity cont = fe->get_continuity();
if (cont == C_ONE)
{
// We'll need gradient data for a C1 projection
libmesh_assert(g.get());
const std::vector<std::vector<RealGradient> >&
ref_dphi = fe->get_dphi();
dphi = &ref_dphi;
}
// The Jacobian * quadrature weight at the quadrature points
const std::vector<Real>& JxW =
fe->get_JxW();
// The XYZ locations of the quadrature points
const std::vector<Point>& xyz_values =
fe->get_xyz();
// The global DOF indices
std::vector<dof_id_type> dof_indices;
// Side/edge DOF indices
std::vector<unsigned int> side_dofs;
// Iterate over all the elements in the range
for (ConstElemRange::const_iterator elem_it=range.begin(); elem_it != range.end(); ++elem_it)
{
const Elem* elem = *elem_it;
// Per-subdomain variables don't need to be projected on
// elements where they're not active
if (!variable.active_on_subdomain(elem->subdomain_id()))
continue;
// Update the DOF indices for this element based on
// the current mesh
dof_map.dof_indices (elem, dof_indices, var);
// The number of DOFs on the element
const unsigned int n_dofs =
cast_int<unsigned int>(dof_indices.size());
// Fixed vs. free DoFs on edge/face projections
std::vector<char> dof_is_fixed(n_dofs, false); // bools
std::vector<int> free_dof(n_dofs, 0);
// The element type
const ElemType elem_type = elem->type();
// The number of nodes on the new element
const unsigned int n_nodes = elem->n_nodes();
// Zero the interpolated values
Ue.resize (n_dofs); Ue.zero();
// In general, we need a series of
// projections to ensure a unique and continuous
// solution. We start by interpolating nodes, then
// hold those fixed and project edges, then
// hold those fixed and project faces, then
// hold those fixed and project interiors
// Interpolate node values first
unsigned int current_dof = 0;
for (unsigned int n=0; n!= n_nodes; ++n)
{
// FIXME: this should go through the DofMap,
// not duplicate dof_indices code badly!
const unsigned int nc =
FEInterface::n_dofs_at_node (dim, fe_type, elem_type,
n);
if (!elem->is_vertex(n))
{
current_dof += nc;
continue;
}
if (cont == DISCONTINUOUS)
{
libmesh_assert_equal_to (nc, 0);
}
// Assume that C_ZERO elements have a single nodal
// value shape function
else if (cont == C_ZERO)
{
libmesh_assert_equal_to (nc, 1);
Ue(current_dof) = f->component(var_component,
elem->point(n),
system.time);
dof_is_fixed[current_dof] = true;
current_dof++;
}
// The hermite element vertex shape functions are weird
else if (fe_type.family == HERMITE)
{
Ue(current_dof) = f->component(var_component,
elem->point(n),
system.time);
dof_is_fixed[current_dof] = true;
current_dof++;
Gradient grad = g->component(var_component,
elem->point(n),
system.time);
// x derivative
Ue(current_dof) = grad(0);
dof_is_fixed[current_dof] = true;
current_dof++;
if (dim > 1)
{
// We'll finite difference mixed derivatives
Point nxminus = elem->point(n),
nxplus = elem->point(n);
nxminus(0) -= TOLERANCE;
nxplus(0) += TOLERANCE;
Gradient gxminus = g->component(var_component,
nxminus,
system.time);
Gradient gxplus = g->component(var_component,
nxplus,
system.time);
// y derivative
Ue(current_dof) = grad(1);
dof_is_fixed[current_dof] = true;
current_dof++;
// xy derivative
Ue(current_dof) = (gxplus(1) - gxminus(1))
/ 2. / TOLERANCE;
dof_is_fixed[current_dof] = true;
current_dof++;
if (dim > 2)
{
// z derivative
Ue(current_dof) = grad(2);
dof_is_fixed[current_dof] = true;
current_dof++;
// xz derivative
Ue(current_dof) = (gxplus(2) - gxminus(2))
/ 2. / TOLERANCE;
dof_is_fixed[current_dof] = true;
current_dof++;
// We need new points for yz
Point nyminus = elem->point(n),
nyplus = elem->point(n);
nyminus(1) -= TOLERANCE;
nyplus(1) += TOLERANCE;
Gradient gyminus = g->component(var_component,
nyminus,
system.time);
Gradient gyplus = g->component(var_component,
nyplus,
system.time);
// xz derivative
Ue(current_dof) = (gyplus(2) - gyminus(2))
/ 2. / TOLERANCE;
dof_is_fixed[current_dof] = true;
current_dof++;
// Getting a 2nd order xyz is more tedious
Point nxmym = elem->point(n),
nxmyp = elem->point(n),
nxpym = elem->point(n),
nxpyp = elem->point(n);
nxmym(0) -= TOLERANCE;
nxmym(1) -= TOLERANCE;
nxmyp(0) -= TOLERANCE;
nxmyp(1) += TOLERANCE;
nxpym(0) += TOLERANCE;
nxpym(1) -= TOLERANCE;
nxpyp(0) += TOLERANCE;
nxpyp(1) += TOLERANCE;
Gradient gxmym = g->component(var_component,
nxmym,
system.time);
Gradient gxmyp = g->component(var_component,
nxmyp,
system.time);
Gradient gxpym = g->component(var_component,
nxpym,
system.time);
Gradient gxpyp = g->component(var_component,
nxpyp,
system.time);
Number gxzplus = (gxpyp(2) - gxmyp(2))
/ 2. / TOLERANCE;
Number gxzminus = (gxpym(2) - gxmym(2))
/ 2. / TOLERANCE;
// xyz derivative
Ue(current_dof) = (gxzplus - gxzminus)
/ 2. / TOLERANCE;
dof_is_fixed[current_dof] = true;
current_dof++;
}
}
}
// Assume that other C_ONE elements have a single nodal
// value shape function and nodal gradient component
// shape functions
else if (cont == C_ONE)
{
libmesh_assert_equal_to (nc, 1 + dim);
Ue(current_dof) = f->component(var_component,
elem->point(n),
system.time);
dof_is_fixed[current_dof] = true;
current_dof++;
Gradient grad = g->component(var_component,
elem->point(n),
system.time);
for (unsigned int i=0; i!= dim; ++i)
{
Ue(current_dof) = grad(i);
dof_is_fixed[current_dof] = true;
current_dof++;
}
}
else
libmesh_error_msg("Unknown continuity " << cont);
}
// In 3D, project any edge values next
if (dim > 2 && cont != DISCONTINUOUS)
for (unsigned int e=0; e != elem->n_edges(); ++e)
{
FEInterface::dofs_on_edge(elem, dim, fe_type, e,
side_dofs);
// Some edge dofs are on nodes and already
// fixed, others are free to calculate
unsigned int free_dofs = 0;
for (unsigned int i=0; i != side_dofs.size(); ++i)
if (!dof_is_fixed[side_dofs[i]])
free_dof[free_dofs++] = i;
// There may be nothing to project
if (!free_dofs)
continue;
Ke.resize (free_dofs, free_dofs); Ke.zero();
Fe.resize (free_dofs); Fe.zero();
// The new edge coefficients
DenseVector<Number> Uedge(free_dofs);
// Initialize FE data on the edge
fe->attach_quadrature_rule (qedgerule.get());
fe->edge_reinit (elem, e);
const unsigned int n_qp = qedgerule->n_points();
// Loop over the quadrature points
for (unsigned int qp=0; qp<n_qp; qp++)
{
// solution at the quadrature point
Number fineval = f->component(var_component,
xyz_values[qp],
system.time);
// solution grad at the quadrature point
Gradient finegrad;
if (cont == C_ONE)
finegrad = g->component(var_component,
xyz_values[qp],
system.time);
// Form edge projection matrix
for (unsigned int sidei=0, freei=0;
sidei != side_dofs.size(); ++sidei)
{
unsigned int i = side_dofs[sidei];
// fixed DoFs aren't test functions
if (dof_is_fixed[i])
continue;
for (unsigned int sidej=0, freej=0;
sidej != side_dofs.size(); ++sidej)
{
unsigned int j = side_dofs[sidej];
if (dof_is_fixed[j])
Fe(freei) -= phi[i][qp] * phi[j][qp] *
JxW[qp] * Ue(j);
else
Ke(freei,freej) += phi[i][qp] *
phi[j][qp] * JxW[qp];
if (cont == C_ONE)
{
if (dof_is_fixed[j])
Fe(freei) -= ((*dphi)[i][qp] *
(*dphi)[j][qp]) *
JxW[qp] * Ue(j);
else
Ke(freei,freej) += ((*dphi)[i][qp] *
(*dphi)[j][qp])
* JxW[qp];
}
if (!dof_is_fixed[j])
freej++;
}
Fe(freei) += phi[i][qp] * fineval * JxW[qp];
if (cont == C_ONE)
Fe(freei) += (finegrad * (*dphi)[i][qp]) *
JxW[qp];
freei++;
}
}
Ke.cholesky_solve(Fe, Uedge);
// Transfer new edge solutions to element
for (unsigned int i=0; i != free_dofs; ++i)
{
Number &ui = Ue(side_dofs[free_dof[i]]);
libmesh_assert(std::abs(ui) < TOLERANCE ||
std::abs(ui - Uedge(i)) < TOLERANCE);
ui = Uedge(i);
dof_is_fixed[side_dofs[free_dof[i]]] = true;
}
}
// Project any side values (edges in 2D, faces in 3D)
if (dim > 1 && cont != DISCONTINUOUS)
for (unsigned int s=0; s != elem->n_sides(); ++s)
{
FEInterface::dofs_on_side(elem, dim, fe_type, s,
side_dofs);
// Some side dofs are on nodes/edges and already
// fixed, others are free to calculate
unsigned int free_dofs = 0;
for (unsigned int i=0; i != side_dofs.size(); ++i)
if (!dof_is_fixed[side_dofs[i]])
free_dof[free_dofs++] = i;
// There may be nothing to project
if (!free_dofs)
continue;
Ke.resize (free_dofs, free_dofs); Ke.zero();
Fe.resize (free_dofs); Fe.zero();
// The new side coefficients
DenseVector<Number> Uside(free_dofs);
// Initialize FE data on the side
fe->attach_quadrature_rule (qsiderule.get());
fe->reinit (elem, s);
const unsigned int n_qp = qsiderule->n_points();
// Loop over the quadrature points
for (unsigned int qp=0; qp<n_qp; qp++)
{
// solution at the quadrature point
Number fineval = f->component(var_component,
xyz_values[qp],
system.time);
// solution grad at the quadrature point
Gradient finegrad;
if (cont == C_ONE)
finegrad = g->component(var_component,
xyz_values[qp],
system.time);
// Form side projection matrix
for (unsigned int sidei=0, freei=0;
sidei != side_dofs.size(); ++sidei)
{
unsigned int i = side_dofs[sidei];
// fixed DoFs aren't test functions
if (dof_is_fixed[i])
continue;
for (unsigned int sidej=0, freej=0;
sidej != side_dofs.size(); ++sidej)
{
unsigned int j = side_dofs[sidej];
if (dof_is_fixed[j])
Fe(freei) -= phi[i][qp] * phi[j][qp] *
JxW[qp] * Ue(j);
else
Ke(freei,freej) += phi[i][qp] *
phi[j][qp] * JxW[qp];
if (cont == C_ONE)
{
if (dof_is_fixed[j])
Fe(freei) -= ((*dphi)[i][qp] *
(*dphi)[j][qp]) *
JxW[qp] * Ue(j);
else
Ke(freei,freej) += ((*dphi)[i][qp] *
(*dphi)[j][qp])
* JxW[qp];
}
if (!dof_is_fixed[j])
freej++;
}
Fe(freei) += (fineval * phi[i][qp]) * JxW[qp];
if (cont == C_ONE)
Fe(freei) += (finegrad * (*dphi)[i][qp]) *
JxW[qp];
freei++;
}
}
Ke.cholesky_solve(Fe, Uside);
// Transfer new side solutions to element
for (unsigned int i=0; i != free_dofs; ++i)
{
Number &ui = Ue(side_dofs[free_dof[i]]);
libmesh_assert(std::abs(ui) < TOLERANCE ||
std::abs(ui - Uside(i)) < TOLERANCE);
ui = Uside(i);
dof_is_fixed[side_dofs[free_dof[i]]] = true;
}
}
// Project the interior values, finally
// Some interior dofs are on nodes/edges/sides and
// already fixed, others are free to calculate
unsigned int free_dofs = 0;
for (unsigned int i=0; i != n_dofs; ++i)
if (!dof_is_fixed[i])
free_dof[free_dofs++] = i;
// There may be nothing to project
if (free_dofs)
{
Ke.resize (free_dofs, free_dofs); Ke.zero();
Fe.resize (free_dofs); Fe.zero();
// The new interior coefficients
DenseVector<Number> Uint(free_dofs);
// Initialize FE data
fe->attach_quadrature_rule (qrule.get());
fe->reinit (elem);
const unsigned int n_qp = qrule->n_points();
// Loop over the quadrature points
for (unsigned int qp=0; qp<n_qp; qp++)
{
// solution at the quadrature point
Number fineval = f->component(var_component,
xyz_values[qp],
system.time);
// solution grad at the quadrature point
Gradient finegrad;
if (cont == C_ONE)
finegrad = g->component(var_component,
xyz_values[qp],
system.time);
// Form interior projection matrix
for (unsigned int i=0, freei=0; i != n_dofs; ++i)
{
// fixed DoFs aren't test functions
if (dof_is_fixed[i])
continue;
for (unsigned int j=0, freej=0; j != n_dofs; ++j)
{
if (dof_is_fixed[j])
Fe(freei) -= phi[i][qp] * phi[j][qp] * JxW[qp]
* Ue(j);
else
Ke(freei,freej) += phi[i][qp] * phi[j][qp] *
JxW[qp];
if (cont == C_ONE)
{
if (dof_is_fixed[j])
Fe(freei) -= ((*dphi)[i][qp] *
(*dphi)[j][qp]) * JxW[qp] *
Ue(j);
else
Ke(freei,freej) += ((*dphi)[i][qp] *
(*dphi)[j][qp]) *
JxW[qp];
}
if (!dof_is_fixed[j])
freej++;
}
Fe(freei) += phi[i][qp] * fineval * JxW[qp];
if (cont == C_ONE)
Fe(freei) += (finegrad * (*dphi)[i][qp]) * JxW[qp];
freei++;
}
}
Ke.cholesky_solve(Fe, Uint);
// Transfer new interior solutions to element
for (unsigned int i=0; i != free_dofs; ++i)
{
Number &ui = Ue(free_dof[i]);
libmesh_assert(std::abs(ui) < TOLERANCE ||
std::abs(ui - Uint(i)) < TOLERANCE);
ui = Uint(i);
dof_is_fixed[free_dof[i]] = true;
}
} // if there are free interior dofs
// Make sure every DoF got reached!
for (unsigned int i=0; i != n_dofs; ++i)
libmesh_assert(dof_is_fixed[i]);
const dof_id_type
first = new_vector.first_local_index(),
last = new_vector.last_local_index();
// Lock the new_vector since it is shared among threads.
{
Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);
for (unsigned int i = 0; i < n_dofs; i++)
// We may be projecting a new zero value onto
// an old nonzero approximation - RHS
// if (Ue(i) != 0.)
if ((dof_indices[i] >= first) &&
(dof_indices[i] < last))
{
new_vector.set(dof_indices[i], Ue(i));
}
}
} // end elem loop
} // end variables loop
}
Member Data Documentation
UniquePtr<FunctionBase<Number> > libMesh::ProjectSolution::f [private] |
Definition at line 115 of file system_projection.C.
Referenced by ProjectSolution().
UniquePtr<FunctionBase<Gradient> > libMesh::ProjectSolution::g [private] |
Definition at line 116 of file system_projection.C.
Referenced by ProjectSolution().
NumericVector<Number>& libMesh::ProjectSolution::new_vector [private] |
Definition at line 118 of file system_projection.C.
const Parameters& libMesh::ProjectSolution::parameters [private] |
Definition at line 117 of file system_projection.C.
const System& libMesh::ProjectSolution::system [private] |
Definition at line 113 of file system_projection.C.
Referenced by operator()().
The documentation for this class was generated from the following file:
generated by
