void distort (MeshBase &mesh, const Real factor, const bool perturb_boundary=false)
void translate (MeshBase &mesh, const Real xt=0., const Real yt=0., const Real zt=0.)
void rotate (MeshBase &mesh, const Real phi, const Real theta=0., const Real psi=0.)
void scale (MeshBase &mesh, const Real xs, const Real ys=0., const Real zs=0.)
void all_tri (MeshBase &mesh)
void smooth (MeshBase &, unsigned int, Real)
void flatten (MeshBase &mesh)
void change_boundary_id (MeshBase &mesh, const short int old_id, const short int new_id)
Tools for Mesh modification.
Author:
Date:
Version:
Revision:
Definition at line 626 of file mesh_modification.C.
References MeshBase::add_elem(), MeshBase::add_point(), MeshBase::boundary_info, MeshBase::delete_elem(), MeshBase::elements_begin(), MeshBase::elements_end(), Utility::enum_to_string< ElemType >(), BoundaryInfo::invalid_id, MeshBase::n_active_elem(), MeshBase::node(), MeshBase::prepare_for_use(), DofObject::processor_id(), libMeshEnums::QUAD4, libMeshEnums::QUAD8, libMeshEnums::QUAD9, Elem::set_node(), Elem::subdomain_id(), libMeshEnums::TRI3, and libMeshEnums::TRI6.
{
// We store pointers to the newly created elements in a vector
// until they are ready to be added to the mesh. This is because
// adding new elements on the fly can cause reallocation and invalidation
// of existing iterators.
std::vector<Elem*> new_elements;
new_elements.reserve (2*mesh.n_active_elem());
// If the original mesh has boundary data, we carry that over
// to the new mesh with triangular elements.
const bool mesh_has_boundary_data = (mesh.boundary_info->n_boundary_ids() > 0);
// Temporary vectors to store the new boundary element pointers, side numbers, and boundary ids
std::vector<Elem*> new_bndry_elements;
std::vector<unsigned short int> new_bndry_sides;
std::vector<short int> new_bndry_ids;
// Iterate over the elements, splitting QUADS into
// pairs of conforming triangles.
{
MeshBase::element_iterator el = mesh.elements_begin();
const MeshBase::element_iterator end = mesh.elements_end();
for (; el!=end; ++el)
{
const ElemType etype = (*el)->type();
// all_tri currently only works on coarse meshes
libmesh_assert ((*el)->parent() == NULL);
// We split the quads using the shorter of the two diagonals
// to maintain the best angle properties.
bool edge_swap = false;
// True if we actually split the current element.
bool split_elem = false;
// The two new triangular elements we will split the quad into.
Elem* tri0 = NULL;
Elem* tri1 = NULL;
switch (etype)
{
case QUAD4:
{
split_elem = true;
tri0 = new Tri3;
tri1 = new Tri3;
// Check for possible edge swap
if (((*el)->point(0) - (*el)->point(2)).size() <
((*el)->point(1) - (*el)->point(3)).size())
{
tri0->set_node(0) = (*el)->get_node(0);
tri0->set_node(1) = (*el)->get_node(1);
tri0->set_node(2) = (*el)->get_node(2);
tri1->set_node(0) = (*el)->get_node(0);
tri1->set_node(1) = (*el)->get_node(2);
tri1->set_node(2) = (*el)->get_node(3);
}
else
{
edge_swap=true;
tri0->set_node(0) = (*el)->get_node(0);
tri0->set_node(1) = (*el)->get_node(1);
tri0->set_node(2) = (*el)->get_node(3);
tri1->set_node(0) = (*el)->get_node(1);
tri1->set_node(1) = (*el)->get_node(2);
tri1->set_node(2) = (*el)->get_node(3);
}
break;
}
case QUAD8:
{
split_elem = true;
tri0 = new Tri6;
tri1 = new Tri6;
Node* new_node = mesh.add_point((mesh.node((*el)->node(0)) +
mesh.node((*el)->node(1)) +
mesh.node((*el)->node(2)) +
mesh.node((*el)->node(3)))*.25
);
// Check for possible edge swap
if (((*el)->point(0) - (*el)->point(2)).size() <
((*el)->point(1) - (*el)->point(3)).size())
{
tri0->set_node(0) = (*el)->get_node(0);
tri0->set_node(1) = (*el)->get_node(1);
tri0->set_node(2) = (*el)->get_node(2);
tri0->set_node(3) = (*el)->get_node(4);
tri0->set_node(4) = (*el)->get_node(5);
tri0->set_node(5) = new_node;
tri1->set_node(0) = (*el)->get_node(0);
tri1->set_node(1) = (*el)->get_node(2);
tri1->set_node(2) = (*el)->get_node(3);
tri1->set_node(3) = new_node;
tri1->set_node(4) = (*el)->get_node(6);
tri1->set_node(5) = (*el)->get_node(7);
}
else
{
edge_swap=true;
tri0->set_node(0) = (*el)->get_node(3);
tri0->set_node(1) = (*el)->get_node(0);
tri0->set_node(2) = (*el)->get_node(1);
tri0->set_node(3) = (*el)->get_node(7);
tri0->set_node(4) = (*el)->get_node(4);
tri0->set_node(5) = new_node;
tri1->set_node(0) = (*el)->get_node(1);
tri1->set_node(1) = (*el)->get_node(2);
tri1->set_node(2) = (*el)->get_node(3);
tri1->set_node(3) = (*el)->get_node(5);
tri1->set_node(4) = (*el)->get_node(6);
tri1->set_node(5) = new_node;
}
break;
}
case QUAD9:
{
split_elem = true;
tri0 = new Tri6;
tri1 = new Tri6;
// Check for possible edge swap
if (((*el)->point(0) - (*el)->point(2)).size() <
((*el)->point(1) - (*el)->point(3)).size())
{
tri0->set_node(0) = (*el)->get_node(0);
tri0->set_node(1) = (*el)->get_node(1);
tri0->set_node(2) = (*el)->get_node(2);
tri0->set_node(3) = (*el)->get_node(4);
tri0->set_node(4) = (*el)->get_node(5);
tri0->set_node(5) = (*el)->get_node(8);
tri1->set_node(0) = (*el)->get_node(0);
tri1->set_node(1) = (*el)->get_node(2);
tri1->set_node(2) = (*el)->get_node(3);
tri1->set_node(3) = (*el)->get_node(8);
tri1->set_node(4) = (*el)->get_node(6);
tri1->set_node(5) = (*el)->get_node(7);
}
else
{
edge_swap=true;
tri0->set_node(0) = (*el)->get_node(0);
tri0->set_node(1) = (*el)->get_node(1);
tri0->set_node(2) = (*el)->get_node(3);
tri0->set_node(3) = (*el)->get_node(4);
tri0->set_node(4) = (*el)->get_node(8);
tri0->set_node(5) = (*el)->get_node(7);
tri1->set_node(0) = (*el)->get_node(1);
tri1->set_node(1) = (*el)->get_node(2);
tri1->set_node(2) = (*el)->get_node(3);
tri1->set_node(3) = (*el)->get_node(5);
tri1->set_node(4) = (*el)->get_node(6);
tri1->set_node(5) = (*el)->get_node(8);
}
break;
}
// No need to split elements that are already triangles
case TRI3:
case TRI6:
continue;
// Try to ignore non-2D elements for now
default:
{
std::cerr << 'Warning, encountered non-2D element '
<< Utility::enum_to_string<ElemType>(etype)
<< ' in MeshTools::Modification::all_tri(), hope that's OK...'
<< std::endl;
}
} // end switch (etype)
if (split_elem)
{
// Be sure the correct ID's are also set for tri0 and
// tri1.
tri0->processor_id() = (*el)->processor_id();
tri0->subdomain_id() = (*el)->subdomain_id();
tri1->processor_id() = (*el)->processor_id();
tri1->subdomain_id() = (*el)->subdomain_id();
if (mesh_has_boundary_data)
{
for (unsigned int sn=0; sn<(*el)->n_sides(); ++sn)
{
short int b_id = mesh.boundary_info->boundary_id(*el, sn);
if (b_id != BoundaryInfo::invalid_id)
{
// Add the boundary ID to the list of new boundary ids
new_bndry_ids.push_back(b_id);
// Convert the boundary side information of the old element to
// boundary side information for the new element.
if (!edge_swap)
{
switch (sn)
{
case 0:
{
// New boundary side is Tri 0, side 0
new_bndry_elements.push_back(tri0);
new_bndry_sides.push_back(0);
break;
}
case 1:
{
// New boundary side is Tri 0, side 1
new_bndry_elements.push_back(tri0);
new_bndry_sides.push_back(1);
break;
}
case 2:
{
// New boundary side is Tri 1, side 1
new_bndry_elements.push_back(tri1);
new_bndry_sides.push_back(1);
break;
}
case 3:
{
// New boundary side is Tri 1, side 2
new_bndry_elements.push_back(tri1);
new_bndry_sides.push_back(2);
break;
}
default:
{
std::cerr << 'Quad4/8/9 cannot have more than 4 sides.' << std::endl;
libmesh_error();
}
}
}
else // edge_swap==true
{
switch (sn)
{
case 0:
{
// New boundary side is Tri 0, side 0
new_bndry_elements.push_back(tri0);
new_bndry_sides.push_back(0);
break;
}
case 1:
{
// New boundary side is Tri 1, side 0
new_bndry_elements.push_back(tri1);
new_bndry_sides.push_back(0);
break;
}
case 2:
{
// New boundary side is Tri 1, side 1
new_bndry_elements.push_back(tri1);
new_bndry_sides.push_back(1);
break;
}
case 3:
{
// New boundary side is Tri 0, side 2
new_bndry_elements.push_back(tri0);
new_bndry_sides.push_back(2);
break;
}
default:
{
std::cerr << 'Quad4/8/9 cannot have more than 4 sides.' << std::endl;
libmesh_error();
}
}
} // end edge_swap==true
} // end if (b_id != BoundaryInfo::invalid_id)
} // end for loop over sides
// Remove the original element from the BoundaryInfo structure.
mesh.boundary_info->remove(*el);
} // end if (mesh_has_boundary_data)
// Add the newly-created triangles to the temporary vector of new elements.
new_elements.push_back(tri0);
new_elements.push_back(tri1);
// Delete the original element
mesh.delete_elem(*el);
} // end if (split_elem)
} // End for loop over elements
}
// Now, iterate over the new elements vector, and add them each to
// the Mesh.
{
std::vector<Elem*>::iterator el = new_elements.begin();
const std::vector<Elem*>::iterator end = new_elements.end();
for (; el != end; ++el)
{
mesh.add_elem(*el);
}
}
if (mesh_has_boundary_data)
{
// By this time, we should have removed all of the original boundary sides
// - except on a hybrid mesh, where we can't 'start from a blank slate'! - RHS
// libmesh_assert (mesh.boundary_info->n_boundary_conds()==0);
// Clear the boundary info, to be sure and start from a blank slate.
// mesh.boundary_info->clear();
// If the old mesh had boundary data, the new mesh better have some.
libmesh_assert (new_bndry_elements.size() > 0);
// We should also be sure that the lengths of the new boundary data vectors
// are all the same.
libmesh_assert (new_bndry_elements.size() == new_bndry_sides.size());
libmesh_assert (new_bndry_sides.size() == new_bndry_ids.size());
// Add the new boundary info to the mesh
for (unsigned int s=0; s<new_bndry_elements.size(); ++s)
mesh.boundary_info->add_side(new_bndry_elements[s],
new_bndry_sides[s],
new_bndry_ids[s]);
}
// Prepare the newly created mesh for use.
mesh.prepare_for_use();
// Let the new_elements and new_bndry_elements vectors go out of scope.
}
Definition at line 1279 of file mesh_modification.C.
References MeshBase::boundary_info, MeshBase::level_elements_begin(), MeshBase::level_elements_end(), and Elem::n_sides().
{
// Only level-0 elements store BCs. Loop over them.
MeshBase::element_iterator el = mesh.level_elements_begin(0);
const MeshBase::element_iterator end_el = mesh.level_elements_end(0);
for (; el != end_el; ++el)
{
Elem *elem = *el;
unsigned int n_sides = elem->n_sides();
for (unsigned int s=0; s != n_sides; ++s)
if (mesh.boundary_info->boundary_id(elem, s) == old_id)
{
mesh.boundary_info->remove_side(elem, s, old_id);
mesh.boundary_info->add_side(elem, s, new_id);
}
}
}
Definition at line 43 of file mesh_modification.C.
References MeshBase::active_elements_begin(), MeshBase::active_elements_end(), MeshTools::find_boundary_nodes(), MeshBase::mesh_dimension(), std::min(), MeshBase::n_elem(), MeshBase::n_nodes(), MeshBase::node(), and TypeVector< T >::unit().
{
libmesh_assert (mesh.n_nodes());
libmesh_assert (mesh.n_elem());
libmesh_assert ((factor >= 0.) && (factor <= 1.));
START_LOG('distort()', 'MeshTools::Modification');
// First find nodes on the boundary and flag them
// so that we don't move them
// on_boundary holds false (not on boundary) and true (on boundary)
std::vector<bool> on_boundary (mesh.n_nodes(), false);
if (!perturb_boundary) MeshTools::find_boundary_nodes (mesh, on_boundary);
// Now calculate the minimum distance to
// neighboring nodes for each node.
// hmin holds these distances.
std::vector<float> hmin (mesh.n_nodes(), 1.e20);
MeshBase::element_iterator el = mesh.active_elements_begin();
const MeshBase::element_iterator end = mesh.active_elements_end();
for (; el!=end; ++el)
for (unsigned int n=0; n<(*el)->n_nodes(); n++)
hmin[(*el)->node(n)] = std::min(hmin[(*el)->node(n)],
static_cast<float>((*el)->hmin()));
// Now actually move the nodes
{
const unsigned int seed = 123456;
// seed the random number generator
std::srand(seed);
// If the node is on the boundary or
// the node is not used by any element (hmin[n]<1.e20)
// then we should not move it.
// [Note: Testing for (in)equality might be wrong
// (different types, namely float and double)]
for (unsigned int n=0; n<mesh.n_nodes(); n++)
if (!on_boundary[n] && (hmin[n] < 1.e20) )
{
// the direction, random but unit normalized
Point dir( static_cast<Real>(std::rand())/static_cast<Real>(RAND_MAX),
(mesh.mesh_dimension() > 1) ?
static_cast<Real>(std::rand())/static_cast<Real>(RAND_MAX)
: 0.,
((mesh.mesh_dimension() == 3) ?
static_cast<Real>(std::rand())/static_cast<Real>(RAND_MAX)
: 0.)
);
dir(0) = (dir(0)-.5)*2.;
if (mesh.mesh_dimension() > 1)
dir(1) = (dir(1)-.5)*2.;
if (mesh.mesh_dimension() == 3)
dir(2) = (dir(2)-.5)*2.;
dir = dir.unit();
mesh.node(n)(0) += dir(0)*factor*hmin[n];
if (mesh.mesh_dimension() > 1)
mesh.node(n)(1) += dir(1)*factor*hmin[n];
if (mesh.mesh_dimension() == 3)
mesh.node(n)(2) += dir(2)*factor*hmin[n];
}
}
// All done
STOP_LOG('distort()', 'MeshTools::Modification');
}
Definition at line 1174 of file mesh_modification.C.
References MeshBase::active_elements_begin(), MeshBase::active_elements_end(), MeshBase::add_elem(), MeshBase::boundary_info, Elem::build(), MeshBase::delete_elem(), MeshBase::elements_begin(), MeshBase::elements_end(), Elem::get_node(), BoundaryInfo::invalid_id, MeshBase::n_active_elem(), Elem::n_nodes(), Elem::n_sides(), Elem::neighbor(), MeshBase::prepare_for_use(), DofObject::processor_id(), Elem::set_node(), Elem::subdomain_id(), and Elem::type().
{
// Algorithm:
// .) For each active element in the mesh: construct a
// copy which is the same in every way *except* it is
// a level 0 element. Store the pointers to these in
// a separate vector. Save any boundary information as well.
// Delete the active element from the mesh.
// .) Loop over all (remaining) elements in the mesh, delete them.
// .) Add the level-0 copies back to the mesh
// Temporary storage for new element pointers
std::vector<Elem*> new_elements;
// BoundaryInfo Storage for element ids, sides, and BC ids
std::vector<Elem*> saved_boundary_elements;
std::vector<short int> saved_bc_ids;
std::vector<unsigned short int> saved_bc_sides;
// Reserve a reasonable amt. of space for each
new_elements.reserve(mesh.n_active_elem());
saved_boundary_elements.reserve(mesh.boundary_info->n_boundary_conds());
saved_bc_ids.reserve(mesh.boundary_info->n_boundary_conds());
saved_bc_sides.reserve(mesh.boundary_info->n_boundary_conds());
{
MeshBase::element_iterator it = mesh.active_elements_begin();
const MeshBase::element_iterator end = mesh.active_elements_end();
for (; it != end; ++it)
{
Elem* elem = *it;
// Make a new element of the same type
Elem* copy = Elem::build(elem->type()).release();
// Set node pointers (they point to nodes in the original mesh)
// The copy's element ID will be set upon adding it to the mesh
for(unsigned int n=0; n<elem->n_nodes(); n++)
copy->set_node(n) = elem->get_node(n);
// Copy over ids
copy->processor_id() = elem->processor_id();
copy->subdomain_id() = elem->subdomain_id();
// This element could have boundary info as well. We need
// to save the (elem, side, bc_id) triples
for (unsigned int s=0; s<elem->n_sides(); s++)
if (elem->neighbor(s) == NULL)
{
short int bc_id = mesh.boundary_info->boundary_id (elem,s);
if (bc_id != BoundaryInfo::invalid_id)
{
saved_boundary_elements.push_back(copy);
saved_bc_ids.push_back(bc_id);
saved_bc_sides.push_back(s);
}
}
// We're done with this element
mesh.delete_elem(elem);
// But save the copy
new_elements.push_back(copy);
}
// Make sure we saved the same number of boundary conditions
// in each vector.
libmesh_assert (saved_boundary_elements.size() == saved_bc_ids.size());
libmesh_assert (saved_bc_ids.size() == saved_bc_sides.size());
}
// Loop again, delete any remaining elements
{
MeshBase::element_iterator it = mesh.elements_begin();
const MeshBase::element_iterator end = mesh.elements_end();
for (; it != end; ++it)
mesh.delete_elem( *it );
}
// Add the copied (now level-0) elements back to the mesh
{
for (std::vector<Elem*>::iterator it = new_elements.begin();
it != new_elements.end();
++it)
mesh.add_elem(*it);
}
// Finally, also add back the saved boundary information
for (unsigned int e=0; e<saved_boundary_elements.size(); ++e)
mesh.boundary_info->add_side(saved_boundary_elements[e],
saved_bc_sides[e],
saved_bc_ids[e]);
// Trim unused and renumber nodes and elements
mesh.prepare_for_use();
}
Definition at line 159 of file mesh_modification.C.
References MeshBase::mesh_dimension(), MeshBase::n_nodes(), MeshBase::node(), and libMesh::pi.
{
libmesh_assert (mesh.mesh_dimension() != 1);
const Real pi = std::acos(-1.);
const Real p = -phi/180.*pi;
const Real t = -theta/180.*pi;
const Real s = -psi/180.*pi;
const Real sp = std::sin(p), cp = std::cos(p);
const Real st = std::sin(t), ct = std::cos(t);
const Real ss = std::sin(s), cs = std::cos(s);
// We follow the convention described at http://mathworld.wolfram.com/EulerAngles.html
// (equations 6-14 give the entries of the composite transformation matrix).
// The rotations are performed sequentially about the z, x, and z axes, in that order.
// A positive angle yields a counter-clockwise rotation about the axis in question.
for (unsigned int n=0; n<mesh.n_nodes(); n++)
{
const Point p = mesh.node(n);
const Real x = p(0);
const Real y = p(1);
const Real z = p(2);
mesh.node(n) = Point(( cp*cs-sp*ct*ss)*x + ( sp*cs+cp*ct*ss)*y + (st*ss)*z,
(-cp*ss-sp*ct*cs)*x + (-sp*ss+cp*ct*cs)*y + (st*cs)*z,
( sp*st)*x + (-cp*st)*y + (ct)*z );
}
}
Definition at line 191 of file mesh_modification.C.
References MeshBase::n_nodes(), MeshBase::node(), and MeshBase::spatial_dimension().
Referenced by DenseVector< T >::operator*=(), and DenseMatrix< T >::operator*=().
{
const Real x_scale = xs;
Real y_scale = ys;
Real z_scale = zs;
if (ys == 0.)
{
libmesh_assert (zs == 0.);
y_scale = z_scale = x_scale;
}
// Scale the x coordinate in all dimensions
for (unsigned int n=0; n<mesh.n_nodes(); n++)
mesh.node(n)(0) *= x_scale;
// Only scale the y coordinate in 2 and 3D
if (mesh.spatial_dimension() > 1)
{
for (unsigned int n=0; n<mesh.n_nodes(); n++)
mesh.node(n)(1) *= y_scale;
// Only scale the z coordinate in 3D
if (mesh.spatial_dimension() == 3)
{
for (unsigned int n=0; n<mesh.n_nodes(); n++)
mesh.node(n)(2) *= z_scale;
}
}
}
Author:
Date:
This implementation assumes every element 'side' has only 2 nodes.
Definition at line 992 of file mesh_modification.C.
References TypeVector< T >::add(), TypeVector< T >::add_scaled(), Elem::build_side(), Elem::child(), Elem::embedding_matrix(), MeshTools::find_boundary_nodes(), Elem::get_node(), DofObject::id(), MeshBase::level_elements_begin(), MeshBase::level_elements_end(), MeshBase::mesh_dimension(), Elem::n_children(), MeshTools::n_levels(), Elem::n_neighbors(), Elem::n_nodes(), MeshBase::n_nodes(), Elem::n_second_order_adjacent_vertices(), Elem::n_vertices(), Elem::neighbor(), MeshBase::node(), Elem::parent(), Elem::point(), Utility::pow(), Elem::second_order_adjacent_vertex(), TypeVector< T >::size(), and MeshTools::weight().
{
libmesh_assert (mesh.mesh_dimension() == 2);
/*
* find the boundary nodes
*/
std::vector<bool> on_boundary;
MeshTools::find_boundary_nodes(mesh, on_boundary);
for (unsigned int iter=0; iter<n_iterations; iter++)
{
/*
* loop over the mesh refinement level
*/
unsigned int n_levels = MeshTools::n_levels(mesh);
for (unsigned int refinement_level=0; refinement_level != n_levels;
refinement_level++)
{
// initialize the storage (have to do it on every level to get empty vectors
std::vector<Point> new_positions;
std::vector<Real> weight;
new_positions.resize(mesh.n_nodes());
weight.resize(mesh.n_nodes());
{
/*
* Loop over the elements to calculate new node positions
*/
MeshBase::element_iterator el = mesh.level_elements_begin(refinement_level);
const MeshBase::element_iterator end = mesh.level_elements_end(refinement_level);
for (; el != end; ++el)
{
/*
* Constant handle for the element
*/
const Elem* elem = *el;
/*
* We relax all nodes on level 0 first
* If the element is refined (level > 0), we interpolate the
* parents nodes with help of the embedding matrix
*/
if (refinement_level == 0)
{
for (unsigned int s=0; s<elem->n_neighbors(); s++)
{
/*
* Only operate on sides which are on the
* boundary or for which the current element's
* id is greater than its neighbor's.
* Sides get only built once.
*/
if ((elem->neighbor(s) != NULL) &&
(elem->id() > elem->neighbor(s)->id()) )
{
AutoPtr<Elem> side(elem->build_side(s));
Node* node0 = side->get_node(0);
Node* node1 = side->get_node(1);
Real node_weight = 1.;
// calculate the weight of the nodes
if (power > 0)
{
Point diff = (*node0)-(*node1);
node_weight = std::pow( diff.size(), power );
}
const unsigned int id0 = node0->id(), id1 = node1->id();
new_positions[id0].add_scaled( *node1, node_weight );
new_positions[id1].add_scaled( *node0, node_weight );
weight[id0] += node_weight;
weight[id1] += node_weight;
}
} // element neighbor loop
}
#ifdef LIBMESH_ENABLE_AMR
else // refinement_level > 0
{
/*
* Find the positions of the hanging nodes of refined elements.
* We do this by calculating their position based on the parent
* (one level less refined) element, and the embedding matrix
*/
const Elem* parent = elem->parent();
/*
* find out which child I am
*/
for (unsigned int c=0; c < parent->n_children(); c++)
{
if (parent->child(c) == elem)
{
/*
*loop over the childs (that is, the current elements) nodes
*/
for (unsigned int nc=0; nc < elem->n_nodes(); nc++)
{
/*
* the new position of the node
*/
Point point;
for (unsigned int n=0; n<parent->n_nodes(); n++)
{
/*
* The value from the embedding matrix
*/
const float em_val = parent->embedding_matrix(c,nc,n);
if (em_val != 0.)
point.add_scaled (parent->point(n), em_val);
}
const unsigned int id = elem->get_node(nc)->id();
new_positions[id] = point;
weight[id] = 1.;
}
} // if parent->child == elem
} // for parent->n_children
} // if element refinement_level
#endif // #ifdef LIBMESH_ENABLE_AMR
} // element loop
/*
* finally reposition the vertex nodes
*/
for (unsigned int nid=0; nid<mesh.n_nodes(); ++nid)
if (!on_boundary[nid] && weight[nid] > 0.)
mesh.node(nid) = new_positions[nid]/weight[nid];
}
{
/*
* Now handle the additional second_order nodes by calculating
* their position based on the vertex postitions
* we do a second loop over the level elements
*/
MeshBase::element_iterator el = mesh.level_elements_begin(refinement_level);
const MeshBase::element_iterator end = mesh.level_elements_end(refinement_level);
for (; el != end; ++el)
{
/*
* Constant handle for the element
*/
const Elem* elem = *el;
const unsigned int son_begin = elem->n_vertices();
const unsigned int son_end = elem->n_nodes();
for (unsigned int n=son_begin; n<son_end; n++)
{
const unsigned int n_adjacent_vertices =
elem->n_second_order_adjacent_vertices(n);
Point point;
for (unsigned int v=0; v<n_adjacent_vertices; v++)
point.add(elem->point( elem->second_order_adjacent_vertex(n,v) ));
const unsigned int id = elem->get_node(n)->id();
mesh.node(id) = point/n_adjacent_vertices;
}
}
}
} // refinement_level loop
} // end iteration
}
Definition at line 125 of file mesh_modification.C.
References MeshBase::n_nodes(), and MeshBase::node().
{
const Point p(xt, yt, zt);
for (unsigned int n=0; n<mesh.n_nodes(); n++)
mesh.node(n) += p;
}
Generated automatically by Doxygen for libMesh from the source code.