Re: [deal.II] Re: mesh quality during deformations using MappingQEulerian

2016-12-01 Thread Jean-Paul Pelteret
Hi Tom,
 

> J-P, you (and the entire dealii community) are a paragon of encouragement, 
>> and I appreciate that very much.  
>>
>
Thank you, your kind words do mean a lot!
 

> here's what I came up with yesterday
>

That's a handy piece of code to keep around. Thanks for sharing it!

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.


Its great! Yeah, a complete remesh is an entirely different story, but if 
you have a look in the documentation of FE_Field_Function, these two lines 

  in 
particular, you will be happy to find that this *might* not be too 
difficult to implement either.
 

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

I keep on forgetting that you're working in codim 1. Yes, 
GridIn::read_abaqus will definitely not cater for that case as it doesn't 
distinguish between dim and spacedim, but it might be an easy fix 

.
 

> I will try VectorTools::interpolate_to_different_mesh 
> 
>  
> when I get to my desk this morning.  I tried the FEFieldFunction 
> 
>  
> yesterday, but 
>
> FEFieldFunction > 
> FEFieldFunction(vector_dof_handler, euler_vector) 
>
> is very unhappy - it is not implemented for DoFHandlerType = 
> DoFHandler (and Mapping), rather it is 
> implemented for dim=spacedim problems only.  
>

I'm not sure if there's any specific reason as to why its only templatised 

 
and instantiated 

 for 
the codim 0 case, but its an enhancement that we could look into making. 
This class makes heavy use of the GridTools::find_cell_around_point 
function, and *perhaps* (pure speculation) it may not be sufficiently 
robust codim 1. Perhaps someone more knowledgable on the inner functioning 
/ use of this class would like to comment on this point?

Cheers,
J-P 

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


Re: [deal.II] Re: mesh quality during deformations using MappingQEulerian

2016-12-01 Thread thomas stephens
> 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
>> 
>> moves
>> 
>> 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 vertex_fe(FE_Q(1), spacedim);
>3   DoFHandler vertex_dof_handler(temp_triangulation);
>4   vertex_dof_handler.distribute_dofs(vertex_fe);
>5   Vector vertex_euler_vector(vertex_dof_handler.n_dofs());
>6
>7   std::vector > original_vertices =
> triangulation.get_vertices();
>8   std::vector vertex_touched(temp_triangulat
> ion.n_vertices(),false);
>9
>   10   typename DoFHandler::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::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>   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.  Ev

Re: [deal.II] Re: mesh quality during deformations using MappingQEulerian

2016-11-30 Thread Jean-Paul Pelteret
Hi Tom,

Its great that you're getting somewhere with this!

I have taken some good advice from Oded (it turns out we work one building 
> over from each other!), 
>

What a happy coincidence :-)
 

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

In theory this shouldn't be a problem. This is sort of a "black-box" step.
 

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

Here's a snippet of one of my cubit scripts that may be of use?

# === Export geometry ===
# - Cubit format
cubit.cmd('save as "' + out_mesh.output_cub_file + '" overwrite')
# - Abaqus format
cubit.cmd('set Abaqus precision 6')
cubit.cmd('export Abaqus "' + out_mesh.output_abaqus_file + '" dimension ' 
+ str(in_geom.dim) + ' overwrite')

2. A function that actually 
>> 
>>  
>> moves 
>> 
>>  
>> 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 vertex_fe(FE_Q(1), spacedim);
>3   DoFHandler vertex_dof_handler(temp_triangulation);
>4   vertex_dof_handler.distribute_dofs(vertex_fe);
>5   Vector vertex_euler_vector(vertex_dof_handler.n_dofs());
>6 
>7   std::vector > original_vertices = 
> triangulation.get_vertices();
>8   std::vector 
> vertex_touched(temp_triangulation.n_vertices(),false);
>9 
>   10   typename DoFHandler::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::vertices_per_cell, right?
 

>   15   if (!vertex_touched[cell->vertex_index(v)])
>   16   {
>   17 for (unsigned int k=0; 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.

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 

 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(FE_Q(n),dim). If I'm 
wrong on this then you could wrap your Euler vector in an FEFieldFunction 

 
and just use apply of the VectorTools::interpolate 

 
functions - this may be slow if the triangulation is huge, but it woul

Re: [deal.II] Re: mesh quality during deformations using MappingQEulerian

2016-11-30 Thread thomas stephens
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
 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
> 
> moves
> 
> 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 vertex_fe(FE_Q(1), spacedim);
   3   DoFHandler vertex_dof_handler(temp_triangulation);
   4   vertex_dof_handler.distribute_dofs(vertex_fe);
   5   Vector vertex_euler_vector(vertex_dof_handler.n_dofs());
   6
   7   std::vector > original_vertices =
triangulation.get_vertices();
   8   std::vector
vertex_touched(temp_triangulation.n_vertices(),false);
   9
  10   typename DoFHandler::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; kvertex_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 vertex_fe and initialize it:
   vertex_fe(FE_Q(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 p

Re: [deal.II] Re: mesh quality during deformations using MappingQEulerian

2016-11-16 Thread Wolfgang Bangerth

On 11/16/2016 03:30 PM, thomas stephens wrote:

but I'm not sure it will redistribute mesh points in tangential
directions.  I think there's a second step to this that moves vertices
around.


The important realization is that there are infinitely many ways of 
parametrizing the same surface. That's easy to see if you just think 
about a curve C that you want to describe as a function (x(s),y(s)) 
mapped from a reference domain s \in [0,1]. There are many such 
functions x(s), y(s) that lead to the same curve. Some move along C 
slower in the beginning and faster at the end, some do it the other way 
around. One particular one moves at constant speed -- that's the one 
that uses the arc length (times a constant) for the parameter s.


The same is true with Eulerian mappings: If you move some of the nodes 
tangentially along the surface, it's still the same surface. (At least 
up to the discretization accuracy.)


Best
 W.

--

Wolfgang Bangerth  email: bange...@colostate.edu
   www: http://www.math.colostate.edu/~bangerth/

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


Re: [deal.II] Re: mesh quality during deformations using MappingQEulerian

2016-11-16 Thread thomas stephens

>
>
> In the case that the mesh lives in a higher dimensional space, you also 
> have to enforce -- either as part of the problem formulation, or as a 
> postprocess step for the output of Mesquite -- that the new nodes still 
> need to lie on the geometry as described before. In other words, nodal 
> points can move *within* the surface, but not perpendicular to it. 
>
> It is my understanding that this is precisely what Bonito et al's approach 
enforces - at least in the setting that I will be working in.  Each 
solution I obtain contains a new displacement (velocity field) *and* a new 
mean curvature of the surface.  The approach I am referring to, called 
geometric consistency, uses the identity *h* = -\laplaceBeltrami *x *to 
obtain points that lie on the new surface (at least, up to supporting the 
same mean curvature *h*).  Since I have the mean curvature of the new 
surface at each time step, I can place points x_{n+1} on the (n+1)-st 
surface according to h_{n+1} = -\lapaceBeltrami x_{n+1}, where h_{n+1} is 
an interpolation of h_n onto some new mesh, and the laplace beltrami 
operator is discretized using a mapping that is interpolated onto the new 
mesh.  This is described in the paper I attached higher up in this thread. 
I'm working on implementing this right now, but I'm not sure it will 
redistribute mesh points in tangential directions.  I think there's a 
second step to this that moves vertices around.  


  

> I have no idea whether that is possible within Mesquite, but in the 
> worst case you can always project back to the previous surface. 
>

One of the bullet points on the Mesquite page 
 is: **Improve surface meshes, 
adapt to surface curvature*, which at least sounds promising - I'm looking 
into that this evening.

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


Re: [deal.II] Re: mesh quality during deformations using MappingQEulerian

2016-11-16 Thread Wolfgang Bangerth

On 11/16/2016 02:28 PM, Jean-Paul Pelteret wrote:

1. Defining the triangulation in its desired format (I think this is
documented within the Mesquite documentation)
2. Translating this Euler vector into something thats relevant to the
FESystem for your problem (this is what you're struggling to
conceptualise, but its really not difficult at all).
3. Dealing with hanging nodes (which if I recall correctly Mesquite does
support).



In the case that the mesh lives in a higher dimensional space, you also 
have to enforce -- either as part of the problem formulation, or as a 
postprocess step for the output of Mesquite -- that the new nodes still 
need to lie on the geometry as described before. In other words, nodal 
points can move *within* the surface, but not perpendicular to it.


I have no idea whether that is possible within Mesquite, but in the 
worst case you can always project back to the previous surface.


Best
 W.

--

Wolfgang Bangerth  email: bange...@colostate.edu
   www: http://www.math.colostate.edu/~bangerth/

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


Re: [deal.II] Re: mesh quality during deformations using MappingQEulerian

2016-11-16 Thread Wolfgang Bangerth


Thomas,


I should not have given Mesquite such short shrift in my reply
yesterday.  With regard to the MappingQEulerian approach, I would
imagine Mesquite would take the mapped points as the mesh points
(basically, what my figures show as the mesh), then smooth those out a
bit, and return new vertex locations.  Then I would need to update the
euler_vector to reflect the new vertex positions - this seems to me like
a difficult step since the euler_vector describes displacements of
degrees of freedom rather than vertices.  My difficulty comes from not
truly understanding how to think about degrees of freedom in that
euler_vector in relation to my naive view of the descriptors of a mesh -
i.e. vertices.


The way you ought to see Eulerian mappings is that the geometry is 
simply described by the graph of a piecewise polynomial function. It 
happens that you represent this piecewise polynomial function by a 
vector-valued finite element field, but this field need not have 
anything to do with the solution of the PDE you're trying to solve and 
may in fact use completely different elements. At least in principle, 
the two also don't need to be defined on the same mesh, though that is 
convenient in practice.


What mesh relaxation would do is replace the Eulerian (geometric) field 
described by one solution vector by another solution vector that leads 
to less distortion of cells. That the solution vector may correspond to 
an interpolating finite element that also has degrees of freedom located 
at vertices is really not all that important. Think simply of the 
Eulerian field as a vector field with arrows from every point of the 
reference domain to a corresponding point in the domain you're trying to 
describe, without thinking of the reference domain being subdivided into 
cells.


I hope this helps a bit.

Best
 W.

--

Wolfgang Bangerth  email: bange...@colostate.edu
   www: http://www.math.colostate.edu/~bangerth/

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


Re: [deal.II] Re: mesh quality during deformations using MappingQEulerian

2016-11-15 Thread Wolfgang Bangerth


Thomas,


My initial mesh is not great - here's where I begin:
|

|GridGenerator::hyper_sphere(triangulation,Point(0,0,0),1.0);

triangulation.set_all_manifold_ids(0);

triangulation.refine_global(initial_global_refinement);
|
|/* interpolate the map from triangulation to desired geometry for the first 
time
   a,b,c are the axes of an ellipsoid */|
|VectorTools::interpolate(bending_dof_handler,
   Initial_map_sphere_to_ellipsoid(a,b,c),
   global_euler_vector);
|
const MappingQEulerian,spacedim> mapping (2, dof_handler,
euler_vector);



doubleInitial_map_sphere_to_ellipsoid::value(constPoint&p,constunsignedintcomponent)
 const
{
  doublenorm_p =p.distance(Point(0,0,0));
  // return the difference between the point on the ellipse and the point on
the original geometry
  returnabc_coeffs[component]*p(component)/norm_p -p(component);
}

|

This gets me an ellipsoid, but the corners of the original hyper_sphere are
just translated onto the ellipsoid by euler_vector.  (The color map shows the
mean curvature of the surface.)






I want to deform this shape in response to a scalar function (illustrated by
the new colormap) on the surface, so I adaptively refine and then let the
evolution go forward (the deformation will reduce a global free energy on the
surface) - here's how that looks:






A few timesteps later the surface (and mesh) evolve to look like this: (the
color map shows the mean curvature again)






This eventually crashes as the quadrilaterals become distorted.


It's not my research field, but there should be literature on these sorts of 
phenomena. I think that ultimately, you need to find a way to relax mesh nodes 
in tangential directions.


Best
 W.

--

Wolfgang Bangerth  email: bange...@colostate.edu
   www: http://www.math.colostate.edu/~bangerth/

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