Dear Bhanu,
By Vertex, I meant the geometric quantity only. As I have used first order
> scalar basis functions i.e FE_Q<2> fe(1) to distribute the degrees of
> freedom by dof_handler.distribute_dofs(fe);. I expect 'n_dofs = n_vertices'
> and also the finite element nodes to be residing directly on the geometric
> vertices. So when I use 'vertices = triangulation.get_vertices();' I
> expect the number of elements(Points) in vertices array to be equal to the
> number of nodes(geometric vertices) in my mesh.
>
This is all correct.
> Also I expected these elements (Points of vertices array) to be
> sorted(indexed) in the increasing order of its global_dof_index (which is
> assigned by 'dof_handler.distribute_dofs(fe)' ).
>
However, this is assumption does not hold in general. Although you've
observed that this appears to hold for FE_Q, a change in the internal
implementation of the FE_Q class (which should be irrelevant to the users
of deal.II) could easily nullify this assertion. Each finite element class
is free to order the degrees of freedom however they want, irrespective of
the ordering of the vertices of the triangulation. In fact you will find
that the FE_Q and FE_DGQ classes, which on the face of it seem somewhat
related, order the DoFs quite differently when used within an FESystem. In
particular, a call to get_unit_support_points()
<https://www.dealii.org/8.4.1/doxygen/deal.II/classFiniteElement.html#aa0169a82a73318973ad38e5785acdb27>
might
render a different ordering of the location of the points at which the DoFs
are supported. I've attached a little program to demonstrate this to you,
and hopefully convince you that making any assumptions of the ordering of
DoFs is not such a great idea.
If you're using finite elements with support points at the vertices then
you can safely use vertex_dof_index
<https://www.dealii.org/8.5.0/doxygen/deal.II/classDoFAccessor.html#a5560151b5407e4851d5c1009c7753764>,
which I mentioned previously. This takes into consideration the internal
ordering of the finite element (or finite element system) which currently
your call to
local_dof_indices[v]
does not.
Regards,
Jean-Paul
--
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.
##
# CMake script for the step-1 tutorial program:
##
# Set the name of the project and target:
SET(TARGET vertex_dof_ordering)
# Declare all source files the target consists of:
SET(TARGET_SRC
${TARGET}.cc
# You can specify additional files here!
)
# Usually, you will not need to modify anything beyond this point...
CMAKE_MINIMUM_REQUIRED(VERSION 2.8.8)
FIND_PACKAGE(deal.II 8.4 QUIET
HINTS ${deal.II_DIR} ${DEAL_II_DIR} ../ ../../ $ENV{DEAL_II_DIR}
)
IF(NOT ${deal.II_FOUND})
MESSAGE(FATAL_ERROR "\n"
"*** Could not locate deal.II. ***\n\n"
"You may want to either pass a flag -DDEAL_II_DIR=/path/to/deal.II to
cmake\n"
"or set an environment variable \"DEAL_II_DIR\" that contains this path."
)
ENDIF()
DEAL_II_INITIALIZE_CACHED_VARIABLES()
PROJECT(${TARGET})
DEAL_II_INVOKE_AUTOPILOT()
/* ---------------------------------------------------------------------
*
* Copyright (C) 1999 - 2015 by the deal.II authors
*
* This file is part of the deal.II library.
*
* The deal.II library is free software; you can use it, redistribute
* it, and/or modify it under the terms of the GNU Lesser General
* Public License as published by the Free Software Foundation; either
* version 2.1 of the License, or (at your option) any later version.
* The full text of the license can be found in the file LICENSE at
* the top level of the deal.II distribution.
*
* ---------------------------------------------------------------------
*/
#include <deal.II/grid/tria.h>
#include <deal.II/grid/tria_accessor.h>
#include <deal.II/grid/tria_iterator.h>
#include <deal.II/grid/grid_generator.h>
#include <deal.II/fe/fe_q.h>
#include <deal.II/fe/fe_dgq.h>
#include <deal.II/fe/fe_system.h>
#include <vector>
#include <iostream>
#include <fstream>
#include <cmath>
using namespace dealii;
int main ()
{
const unsigned int dim = 2;
Triangulation<dim> triangulation;
GridGenerator::hyper_cube (triangulation);
const FE_Q<dim> fe_q (1);
const FE_DGQ<dim> fe_dgq (1);
const FESystem<dim> fe_q_sys (fe_q, 2);
const FESystem<dim> fe_dgq_sys (fe_dgq, 2);
for (typename Triangulation<dim>::active_cell_iterator
cell = triangulation.begin_active(); cell != triangulation.end(); ++cell)
{
for (unsigned int v=0; v<GeometryInfo<dim>::vertices_per_cell; ++v)
{
std::cout
<< "TRIA vertex " << v << ": "
<< cell->vertex(v)
<< std::endl;
}
}
for (unsigned int i=0; i<fe_q.get_unit_support_points().size(); ++i)
{
std::cout
<< "FE_Q SP " << i << ": "
<< fe_q.get_unit_support_points()[i]
<< std::endl;
}
for (unsigned int i=0; i<fe_dgq.get_unit_support_points().size(); ++i)
{
std::cout
<< "FE_DGQ SP " << i << ": "
<< fe_dgq.get_unit_support_points()[i]
<< std::endl;
}
for (unsigned int i=0; i<fe_q_sys.get_unit_support_points().size(); ++i)
{
std::cout
<< "FE_Q_SYSTEM SP " << i << ": "
<< fe_q_sys.get_unit_support_points()[i]
<< std::endl;
}
for (unsigned int i=0; i<fe_dgq_sys.get_unit_support_points().size(); ++i)
{
std::cout
<< "FE_DGQ SYSTEM SP " << i << ": "
<< fe_dgq_sys.get_unit_support_points()[i]
<< std::endl;
}
// triangulation.refine_global (4);
// std::ofstream out ("grid-1.eps");
// GridOut grid_out;
// grid_out.write_eps (triangulation, out);
// std::cout << "Grid written to grid-1.eps" << std::endl;
}
