> Its great that you're getting somewhere with this!
>
> J-P, you (and the entire dealii community) are a paragon of encouragement,
and I appreciate that very much.

After I wrote you yesterday I was able to figure out how to 'play a journal
file' from within an embedded python session within dealii, so my problem

  (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)
>

is no longer a problem.  Thank you for your code snippet, and here's what I
came up with yesterday: (this is from within my dealii code)


  // taken from: http://stackoverflow.com/questions/12142174/run-a-
python-script-with-arguments

  char export_script[50];
  char export_filename[150];
  sprintf(export_script,   "export_codim_1_to_ucd.py");
  sprintf(export_filename, "/home/tds3/com/biomembranes/
scratch/mesh_modification/cubit_strategy/data/smoothed_from_cubit.ucd");

  FILE *file;
  file = fopen(export_script, "r");
  size_t argc = 2;
  char *argv[argc];
  argv[0] = export_script;
  argv[1] = export_filename;

  // all this junk requires Python version < 3.0
  // because PySys_SetArgv() does not like the char* to wchar_t* conversion
of argv
  Py_SetProgramName(argv[0]);
  Py_Initialize();

  PyRun_SimpleString("import cubit, sys\n"
                     "cubit.init('')\n"
                     "print('cubit imported')\n"
                     "cubit.cmd(\"import abaqus mesh geometry
'./data/abaqus_test.inp'\")\n"
                     "cubit.cmd(\"smooth surface all laplace\")\n"
                    );
  PySys_SetArgv(argc,argv);
  PyRun_SimpleFile(file, export_script);

  Py_Finalize();

(As a quick aside - I call this export script because I can't get cubit to
export codimension 1 surfaces using its export commands.  The output is
always flattened in the sense that the z-coordinates of vertices are
stripped.  Perhaps your code snippet will help me by setting the dimension
to 3.  For now this works nicely, and cubit is slowly responding to my
email.
As well, I choose ucd because GridIn::read_abaqus() doesn't like to find
that z-vertex from standard abaqus .inp files for codimension 1 surfaces,
but read_ucd() works as expected.)


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.
>


> Ok, so this might be a point of difficulty. See the answer to your next
> question...
>

>
>> 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_triangulat
> ion.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)
>
>
By "4", you mean GeometryInfo<spacedim>::vertices_per_cell, right?
>

YES! Sloppy coding on my part


>   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
>>
>>
>
At first glance this looks ok. My biggest concern is that you've assumed
> that your old and new triangulations have the same vertex ordering. I'm
> inclined to think that this is ok because I don't think that Cubit
> renumbers mesh entities unless you ask it to, and nor does deal.II. But you
> should probably hand-check this on a small geometry.
>

This works on a small geometry, thankfully.  Even after mesh modifications
from within cubit.  This would not work for remeshing from within cubit,
unfortunately.  The complete remeshing case is interesting to me, but
that's not what I'm aimed at right now.

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!
>>
>
I think that VectorTools::interpolate_to_different_mesh
> <https://www.dealii.org/developer/doxygen/deal.II/namespaceVectorTools.html#af68148d58c8dfd0916eceab9d89d74d5>
>  could
> do the job for you; even though it needs the same FE discretisation I
> *think* this means that both DoFHandlers (both based on the original
> triangulation) must be structured like FESystem<dim>(FE_Q(n),dim). If I'm
> wrong on this then you could wrap your Euler vector in an FEFieldFunction
> <https://www.dealii.org/developer/doxygen/deal.II/classFunctions_1_1FEFieldFunction.html>
> and just use apply of the VectorTools::interpolate
> <https://www.dealii.org/developer/doxygen/deal.II/namespaceVectorTools.html#a05db6c8cebf924b417dd92f525efe3db>
> functions - this may be slow if the triangulation is huge, but it would be
> a good first step to take.
>

I will try VectorTools::interpolate_to_different_mesh
<https://www.dealii.org/developer/doxygen/deal.II/namespaceVectorTools.html#af68148d58c8dfd0916eceab9d89d74d5>
when I get to my desk this morning.  I tried the FEFieldFunction
<https://www.dealii.org/developer/doxygen/deal.II/classFunctions_1_1FEFieldFunction.html>
yesterday, but

FEFieldFunction<spacedim, DoFHandler<dim,spacedim> >
FEFieldFunction(vector_dof_handler, euler_vector)

is very unhappy - it is not implemented for DoFHandlerType =
DoFHandler<dim,spacedim> (and Mapping<dim,spacedim>), rather it is
implemented for dim=spacedim problems only.  It's not immediately clear to
me why it is only implemented for dim=spacedim, but it is also not
immediately clear to me how to make that implementation change - perhaps I
will just try to make the obvious changes to that function and see what
happens.


Does any of this sound reasonable?
>>
>
Well, none of it sounds crazy :-)
>

Well, it *feels* a little crazy!!

Thanks so much for your help J-P, and Oded, and Wolfgang (if any of you are
still reading this anymore).  We're almost there, and all of this code is
yours, so hopefully someone else will find it useful.

Tom

-- 
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 dealii+unsubscr...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.
#    Python script to be run as a 'Custom Tool' from the Trelis/Cubit GUI interface. 
#    Exports 3d vertex coordinates and quadrilateral face connectivity
#    information from a Trelis/Cubit mesh on codimension 1 surfaces in
#    what dealii (dealii.org) considers to be the .ucd format for codimension 1 surfaces:
#
#    <number of verts> <number of faces> 0 0 0
#    1 x1 y1 z1
#    2 x2 y2 z2
#    .
#    .
#    .
#    1 <material id> quad q1_v1 q1_v2 q1_v3 q1_v4
#    2 <material id> quad q2_v1 q2_v2 q2_v3 q2_v4 
#    .
#    .
#    .
#
#    where q1_v1 is the first of four vertices associated with quad 1.
#    For example, the codimension 1 surface described by the union of
#    faces of the unit cube
#
#                   4_____________________________________ 3
#                  /|                                    /|
#                 / |             (top) q1              / |
#                /  |                                  /  |
#               /   |                                 /   |
#              /    |                                /    |
#             /     |                               /     |
#           1/____________________________________2/      |
#           |       |                             |       |
#           |       |                             |       |
#           | (side)|                             | (side)|
#           |   q3  |                             |   q5  |
#           |       |                             |       |
#           |       |             (back) q4       |       |
#           |       7-----------------------------|-------8  
#           |      /                              |     / 
#           |     /                 (bottom) q2   |    /  
#           |    /                                |   /   
#           |   /                                 |  /      
#           |  /                                  | /      
#           | /    (front) q6                     |/       
#           6------------------------------------- 5
#
#     has output:
#
#    8 6 0 0 0
#    1 1 0 1
#    2 1 1 1
#    3 0 1 1
#    4 0 0 1
#    5 1 1 0
#    6 1 0 0 
#    7 0 0 0
#    8 0 1 0
#    1 1 quad 1 2 3 4
#    2 1 quad 5 6 7 8
#    3 1 quad 4 7 6 1
#    4 1 quad 3 8 7 4
#    5 1 quad 2 5 8 3
#    6 1 quad 1 6 5 2
#
#    where the <material_id> parameter is the same on all faces, and has been set (abitrarily) to 1



material_id = 1

#output_filename = "/home/tds3/gpde/code/dealii/sandbox/load_mesh/ellipsoid_mesh.ucd"
#output_filename = "/home/tds3/gpde/code/dealii/sandbox/load_mesh/sphere_mesh.ucd"
#output_filename = "/home/tds3/com/biomembranes/scratch/mesh_modification/cubit_strategy/data/smoothed_from_cubit.ucd"
output_filename = sys.argv[0]

#"/home/tds3/com/biomembranes/scratch/mesh_modification/cubit_strategy/data/smoothed_from_cubit.ucd"

# get geometric entities: nodes(vertices) and quads(faces)
nodes = cubit.get_entities("node")
len_nodes = len(nodes)

surfaces = cubit.get_entities("surface")
quads = []
for surface_id in surfaces:
    surface_quads = cubit.get_surface_quads(surface_id)
    quads.extend(surface_quads) 
len_quads = len(quads)

print("nodes: %d, quads %d" %(len_nodes,len_quads))

output_file = open(output_filename,"w")
header = "%d %d 0 0 0\n" %(len_nodes,len_quads)
output_file.write(header)

node_mapping = {} 
for i,node_id in enumerate(nodes):
    coords = cubit.get_nodal_coordinates(node_id)
    coords_str = "%d %0.12f %0.12f %0.12f\n" %(i+1,coords[0],
                                                   coords[1],
                                                   coords[2])
    # map Cubit's node_id's to a more 'natural' enumeration
    node_mapping[node_id] = i+1
    #print(coords_str)
    output_file.write(coords_str)


for j,quad_id in enumerate(quads):
    cubit_connectivity = cubit.get_connectivity("Quad",quad_id)
    # connectivity information comes in the form of Cubit's node_id's,
    # but we want to point back to the enumeration above
    connectivity = (node_mapping[cubit_connectivity[0]],
                    node_mapping[cubit_connectivity[1]], 
                    node_mapping[cubit_connectivity[2]], 
                    node_mapping[cubit_connectivity[3]])
    conn_str = "%d %d quad %d %d %d %d\n" %(j+1,material_id, connectivity[0],
                                                             connectivity[1],
                                                             connectivity[2],
                                                             connectivity[3])
    #print(conn_str)
    output_file.write(conn_str)



print("written %s" % (output_filename))
output_file.close()
print("----- done ---------")

Reply via email to