Dear all,
I would like your guidance on how to perform the assembly of different
shape representations on the same triangulation and on the same loop
through cells.
Let me try to explain a bit more.
I have designed a grid containing two concentric squares (or circles).
Additionally, I have a parametric representation T(alpha) that maps points
(x,y) on the edges
of the inner square (or circle) to points (x',y') on the edges of a desired
shape
contained inside the outer square (or circle). This is
(x',y')=T(alpha;x,y), where
the representation T depends on a parameter alpha.
Now assume we have N different shapes that are labelled with the index m.
We then have N parameters alpha_m, with m=1,2,...N.
Then in the loop through the triangulation cells in the assembly process, I
would like
to be able for each cell to loop through different alpha_m, in order to
generate the
local FE matrices and load vectors corresponding to each of the
modified triangulations corresponding to each of the N shapes.
My guess is that somehow I should modify the following line:
const FEValues<dim> &fe_v = hp_fe_v.get_present_fe_values ();
and pass instead a FEValues that has been generated with the desired
alpha_m.
Below I provide a more complete sketch on how my code looks like.
I am quite hesitant on where to start, and I would appreciate your insights
on
how to achieve having an assemly routine where I can pass different alpha_m
per each
cell in the loop through cells.
I am grateful for any advise or hint on how to achieve this.
Thanks in advance,
Juan Carlos Araújo, PhD
The way I work with one shape is the following:
// Deaclare my environment
PolarShapeManifold manifold;
PolarManifold<dim> polar;
TransfiniteInterpolationManifold<dim> inner_manifold;
Triangulation<dim> triangulation;
hp::DoFHandler<dim> dof_handler;
hp::FECollection<dim> fe_collection;
hp::MappingCollection<dim> mapping_collection;
PETScWrappers::SparseMatrix system_matrix;
PETScWrappers::MPI::Vector solution;
// Constructor
manifold ( alpha), dof_handler (triangulation),...
// setup_system
typename hp::DoFHandler<dim>::active_cell_iterator
cell = dof_handler.begin_active(),
endc = dof_handler.end();
for (; cell!=endc; ++cell) {
cell->set_active_fe_index ( ... );
}
dof_handler.distribute_dofs (fe_collection);
...
// create_coarse_grid
concentric_disks_inner_shape ( triangulation );
tria.set_manifold ( manifold_label, manifold );
tria.set_manifold ( layer_label, polar );
unsigned int not_curved = 1;
inner_manifold.initialize(triangulation);
tria.set_manifold (not_curved, inner_manifold);
...
// Assembly
hp::FEValues<dim> hp_fe_v ( mapping_collection,
fe_collection,
quadrature_collection,
update_values | update_gradients |
update_quadrature_points |
update_JxW_values);
FullMatrix<PetscScalar> cell_system;
Vector<PetscScalar> cell_rhs;
std::vector<types::global_dof_index> local_dof_indices;
typename hp::DoFHandler<dim>::active_cell_iterator
cell = dof_handler.begin_active(),
endc = dof_handler.end();
for (; cell!=endc; ++cell) {
const unsigned int dofs_per_cell = cell->get_fe().dofs_per_cell;
cell_system.reinit (dofs_per_cell, dofs_per_cell);
cell_rhs.reinit (dofs_per_cell);
...
hp_fe_v.reinit (cell);
const FEValues<dim> &fe_v = hp_fe_v.get_present_fe_values ();
for (unsigned int q_point=0; q_point<fe_v.
n_quadrature_points; ++q_point) {
Point<dim> x = fe_v.quadrature_point (q_point);
for (unsigned int i=0; i<dofs_per_cell; ++i){
for (unsigned int j=0; j<dofs_per_cell; ++j) {
cell_system (i,j) += (
fe_v.
shape_grad(i,q_point) *
fe_v.shape_grad(j,q_point)
-f(x) * // given f
fe_v.shape_value(i,q_point) *
fe_v.shape_value(j,q_point) ) *
fe_v.JxW(q_point);
}
cell_rhs (i) += g(x) * // given g
fe_v.
shape_value(i,q_point) *
fe_v.JxW(q_point);
}
} // q
--
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].
To view this discussion on the web visit
https://groups.google.com/d/msgid/dealii/92fc053d-ca4e-4830-bca0-2f211a9ba384%40googlegroups.com.
