[Libmesh-users] Problem with user-generated mesh

[Brent K. <bkr...@ic...>]

Hey all,

I believe there is a bug somewhere in the read-in portion for .xda 
files. I have not been able to track it down, but the code having 
problems. To be sure this wasn't a fixed problem I downloaded the most 
recent version on the sourceforge svn server today.

I generated a tetrahedral mesh using another program and converted it to 
.xda format using a script, based on the answers benkirk gave to my 
question about the format. I have attached the file, but here is the 
header and the first few lines of the connectivity array:


\begin{verbatim}
libMesh-0.7.0+ 
 1143     # number of elements
  222     # number of nodes
.         # boundary condition specification file
n/a       # subdomain id specification file
n/a       # processor id specification file
n/a       # p-level specification file
 1143     # n_elem at level 0, [type (n0 ... nN-1) ]
8     85    19    80    93
8      2     1    26    49
8     77    55    58    17
8     44    57    49    15
\end{verbatim}


I then used the file in a modified version of example 3 (single 
processor) which worked just fine in 3d using using

\begin{verbatim}
MeshTools::Generation::build_cube (mesh,
                                20, 20, 20,
                                       -1., 1.,
                                       -1., 1.,
                                       -1., 1.,
                                       HEX20);
\end{verbatim}


However, when attempting to use this mesh I get the following error 
messages:



\begin{verbatim}
ERROR: Finite element LAGRANGE on geometric element TET
only supports FEInterface::max_order = 1, not fe_type.order = 2
WARNING: Finite element LAGRANGE on geometric element TET
could not be p refined past FEInterface::max_order = 1
ERROR: unsigned char too small to hold Elem::_p_level!
Recompile with Elem:_p_level set to something bigger!
[0] /home/bkraczek/codes-lnx/libmesh/include/geom/elem.h, line 1574, 
compiled Jan 29 2009 at 16:23:47
terminate called after throwing an instance of 'libMesh::LogicError'
  what():  Error in libMesh internal logic
Aborted
\end{verbatim}


What I have found through digging is that when it goes to reinitialize 
the mesh (DofMap::reinit) it finds that fe_type.order == 2;  
elem->p_level()==255 .



Any help would be appreciated,
Brent






Kirk, Benjamin (JSC-EG) wrote:
>> Connectivity:
>> The header is followed immediately by the connectivity. Are blank or
>> commented lines ignored? (i.e. can I put a blank line anywhere?)
>>     
>
> I think so, I'll check...  if not I will fix the code because that is a
> perfectly reasonable request.  Certainly a comment can follow data on a
> line, see for example reference_elements/3D/one_hex20.xda
>
>   
>> In connectivity section, it appears that the first number is the element
>> type. I want to use tetrahedra, so it looks like this number should be 8
>> for a tet4 and 9 for a tet10. The rest of the numbers are the node index
>> numbers, starting from zero.
>>     
>
> Right.
>
>   
>> Nodes:
>> Next come the node locations in Cartesian coordinates. Their order
>> should correspond to the indices used by the connectivity and the
>> possibly the boundary conditions.
>>     
>
> Right.
>
>  
>   
>> Boundary conditions (where I have the most questions):
>>
>> First, is this section necessary if I intend to use the same boundary
>> condition on all external elements? In the examples I've been working
>> from (primarily #3, which uses a penalty method), the code checks to see
>> if the element has a neighbor and applies the BC if it does not. In my
>> application I need to use a special mixed ("Robin") BC, so I don't to
>> use a pure Dirichlet condition.
>>     
>
> Correct.  The boundary conditions are only there for your convenience.  If
> you do not need to 
>
>   
>> Second, how is the side index determined? Specifically, If I have a
>> tetrahedron with arbitrary position in space, how do I know which index
>> a side gets?
>>     
>
> The side indices is dependent on the node ordering.  see for example
> src/geom/cell_tet4.C:
>
> const unsigned int Tet4::side_nodes_map[4][3] =
> {
>   {0, 2, 1}, // Side 0
>   {0, 1, 3}, // Side 1
>   {1, 2, 3}, // Side 2
>   {2, 0, 3}  // Side 3
> };
>
> So side 0 is the triangle defined by nodes (0,2,1) from the tet.
>
>  
>   
>> Third, how is the number for the boundary condition used? Is this index
>> the same as used in something like 'elem->neighbor(side)'?
>>     
>
> It is not used at all.  It is inserted into the BoundaryInfo object
> ('mesh.boundary_info') and you can get access to it to set different BCs on
> different parts of the domain.
>
> -Ben
>
>
>