From: KIRK, BENJAMIN (JSCEG) (NASA) <benjamin.kirk1@na...>  20050728 14:25:48

Please! Ben Original Message From: Martin L=FCthi [mailto:luthi@...]=20 Sent: Wednesday, July 27, 2005 7:10 PM To: KIRK, BENJAMIN (JSCEG) (NASA) Cc: 'Roy Stogner'; Michael Povolotskyi; = libmeshusers@... Subject: RE: [Libmeshusers] how to change node position Michael Some time ago I commited a patch to mesh::smooth, that apparently never = got applied. With this you should be able to move *all* nodes with the = method Ben mentioned: KIRK, BENJAMIN (JSCEG) (NASA) writes: > Functionally, to move a node you simply assign a new Point value to = it: >=20 > mesh.node(n) =3D Point(x_new, y_new, z_new); >=20 > for example. A number of functions in src/mesh/mesh_modification.C = alter > the mesh in such a fashion. After that you make a call to mesh.smooth. What it does is: o relaxes all nodes (you could comment out that step if you don't want this. I used it for an advected mesh, and it helped improve the quality. =20 o looks for higher order nodes (the ones on the element sides) and places them on the edge of the element. This works in 2D, but I'm not sure anymore about 3D. Besides the fact = that the algorithm works for higher order elements it also is much more = efficient than the LaplaceSmoother in the library, since there is no need to = build an expensive data structure of the graph. Hopefully this code still works! Best, Martin BEN: should I prepare another patch? =3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D from = mesh_modification.C = =3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D void MeshTools::Modification::smooth (MeshBase& mesh, const unsigned int n_iterations, const Real power) { /* * find the boundary nodes */ std::vector<bool> on_boundary; MeshTools::find_boundary_nodes(mesh, on_boundary); for (unsigned int iter=3D0; iter<n_iterations; iter++) { /* * loop over the mesh refinement level * TODO: the maximum should be mesh.max_refinement_level(), * but this does not exist (yet) */ for (unsigned int refinement_level=3D0; refinement_level<10; 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 =3D mesh.level_elements_begin(refinement_level); const MeshBase::element_iterator end =3D mesh.level_elements_end(refinement_level);=20 =09 for (; el !=3D end; ++el) { /* * Constant handle for the element */ const Elem* elem =3D *el; =20 /* * We relax all nodes on level 0 first=20 * If the element is refined (level > 0), we interpolate = the * parents nodes with help of the embedding matrix */ if (refinement_level =3D=3D 0) { for (unsigned int s=3D0; 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) !=3D NULL) && (elem>id() > elem>neighbor(s)>id()) )=20 { AutoPtr<Elem> side(elem>build_side(s)); Node* node0 =3D side>get_node(0); Node* node1 =3D side>get_node(1); Real node_weight =3D 1.; // calculate the weight of the nodes if (power > 0) { Point diff =3D (*node0)(*node1); node_weight =3D pow( diff.size(), power = ); } const unsigned int id0 =3D node0>id(), id1 = =3D node1>id(); new_positions[id0].add_scaled( *node1, = node_weight ); new_positions[id1].add_scaled( *node0, = node_weight ); weight[id0] +=3D node_weight; weight[id1] +=3D node_weight; } } // element neighbor loop }=20 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 =3D elem>parent(); /* * find out which child I am */ for (unsigned int c=3D0; c < parent>n_children(); = c++) { if (parent>child(c) =3D=3D elem) { /* *loop over the childs (that is, the current elements) nodes */ for (unsigned int nc=3D0; nc < = elem>n_nodes(); nc++) { /* * the new position of the node */ Point point; for (unsigned int n=3D0; = n<parent>n_nodes(); n++) { /* * The value from the embedding = matrix */ const float em_val =3D parent>embedding_matrix(c,nc,n); =20 if (em_val !=3D 0.) point.add_scaled (parent>point(n), em_val); } const unsigned int id =3D elem>get_node(nc)>id(); new_positions[id] =3D point; weight[id] =3D 1.; } =20 } // if parent>child =3D=3D elem } // for parent>n_children } // if element refinement_level /* * Now handle the additional second_order nodes by = calculating * their position based on the vertex postitions */ std::vector<unsigned int> adjacent_vertices_ids; const unsigned int son_begin =3D elem>n_vertices(); const unsigned int son_end =3D elem>n_nodes(); for (unsigned int n=3Dson_begin; n<son_end; n++) { const unsigned int n_adjacent_vertices =3D elem>n_second_order_adjacent_vertices(n); Point point;=20 for (unsigned int v=3D0; v<n_adjacent_vertices; v++) point.add(elem>point( elem>second_order_adjacent_vertex(n,v) )); const unsigned int id =3D elem>get_node(n)>id(); new_positions[id] =3D point; weight[id] =3D 1.*n_adjacent_vertices; } } // element loop /* * finally reposition the nodes */ for (unsigned int nid=3D0; nid<mesh.n_nodes(); ++nid) if (!on_boundary[nid] && weight[nid] > 0.) mesh.node(nid) =3D new_positions[nid]/weight[nid]; =20 } // refinement_level loop } // end iteration } 