J-P, I have taken some good advice from Oded (it turns out we work one building over from each other!), and I started using Cubit to smooth my meshes. Cubit requires a license but it does seem to have some usage within the dealii community, so I hope whatever comes out of this is useful for somebody. (I have noted your python scripts <https://github.com/dealii/dealii/wiki/Mesh-Input-And-Output> interfacing with Cubit on the githib wiki.)
I'm trying hard to maintain the MappingQEulerian approach in my current model setup, and I wanted to return to the outline you put forward to see if 'Mesquite' can simply be replaced with 'Cubit' in your writeup: My thinking is that one could add the following functionality to the > library: > 1. A function that takes in a triangulation and a (possibly empty) Euler > vector defining the initial mapped coordinates of the triangulation > vertices, and some boolean constraints vector. This would then return an > optimal Q1 Euler vector as computed by Mesquite. > At the moment I have a function that takes an euler_vector defining initial mapped coordinates, along with a dof_handler, for an existing surface and returns a new Triangulation object containing the smoothed mapped vertices from Cubit. The euler_vector is of size dof_handler.n_dofs(), and I would like this to remain arbitrary so I can use higher-order elements (ultimately I'm computing mean curvature of the surface). The Triangulation holding the smoothed vertices arises because I am using GridIn::read_ucd(std::ifsream in) to get data back from cubit. The strategy here is to write my original surface to an abaqus file (code to accomplish this is attached), manipulate it and export it back to dealii using the Python API for Cubit during runtime, and read it back using GridIn. (exporting from cubit is actually a bit of a problem, I have submitted a service email to them - but if I do this through the gui it all works) > 2. A function that actually > <https://www.dealii.org/developer/doxygen/deal.II/step_18.html#TopLevelmove_mesh> > moves > <https://www.dealii.org/developer/doxygen/deal.II/step_42.html#PlasticityContactProblemmove_mesh> > the triangulation vertices for you based on this vector. > At this point the Triangulation attached to my GridIn object is the new triangulation, so this may not be a problem. > 3. A function that would interpolate the optimal Q1 Euler vector onto a > given FE space, which would presumably represent the displacement solution > to some other problem. > > *Here's where my question for you comes in:* At this stage I have a new Triangulation (call it temp_triangulation) which holds my optimal vertices. The vertices in temp_triangulation live on some surface in real space. I have my original Triangulation (call it orig_triangulation) which satisfies orig_triangulation.n_vertices() == temp_triangulation.n_vertices(). *How to I obtain a new euler_vector of size dof_handler.n_dofs() at this stage? * I have done the following so far: 1 /* update euler_vector to describe smoothed mesh */ 2 FESystem<dim,spacedim> vertex_fe(FE_Q<dim,spacedim>(1), spacedim); 3 DoFHandler<dim,spacedim> vertex_dof_handler(temp_triangulation); 4 vertex_dof_handler.distribute_dofs(vertex_fe); 5 Vector<double> vertex_euler_vector(vertex_dof_handler.n_dofs()); 6 7 std::vector<Point<spacedim> > original_vertices = triangulation.get_vertices(); 8 std::vector<bool> vertex_touched(temp_triangulation.n_vertices(),false); 9 10 typename DoFHandler<dim,spacedim>::active_cell_iterator 11 cell=vertex_dof_handler.begin_active(), 12 endc=vertex_dof_handler.end(); 13 for (; cell!=endc; ++cell) 14 for (unsigned int v=0; v < 4; ++v) 15 if (!vertex_touched[cell->vertex_index(v)]) 16 { 17 for (unsigned int k=0; k<spacedim; ++k) 18 { 19 vertex_euler_vector(cell->vertex_dof_index(v,k)) 20 = cell->vertex(v)(k) - original_vertices[cell->vertex_index(v)](k); 21 } 22 vertex_touched[cell->vertex_index(v)] = true; 23 } Basically, I have built vertex_euler_vector to map the vertices of my original reference domain to new points which are more or less on my original surface. This holds essentially Q1 information - only information at the nodes of the original reference domain. I *think* that what I need to do now is interpolate vertex_euler_vector onto my original finite element space (which was Q2, but could be anything). If that is what is needed, I would appreciate some advice on how to accomplish it! Does any of this sound reasonable? - Create an FESystem<dim,spacedim> vertex_fe and initialize it: vertex_fe(FE_Q<dim,spacedim>(1),spacedim) and create a DoFHandler vertex_dof_handler and initialize it with temp_triangulation: vertex_dof_handler(temp_triangulation), and call vertex_dof_handler.distribute_dofs(vertex_fe) - Build a Vector called vertex_euler_vector and reinit it with vertex_dof_handler.n_dofs(): vertex_euler_vector.reinit( vertex_dof_hander.n_dofs()) - Loop through the cells of temp_triangulation and put vertex_euler_vector(vertex_dof_index(v)) = v - orig_triangulation() -- The deal.II project is located at http://www.dealii.org/ For mailing list/forum options, see https://groups.google.com/d/forum/dealii?hl=en --- You received this message because you are subscribed to the Google Groups "deal.II User Group" group. To unsubscribe from this group and stop receiving emails from it, send an email to [email protected]. For more options, visit https://groups.google.com/d/optout.
