From: Wolfgang Bangerth Date: Wed, 26 Apr 2006 04:49:49 +0000 (+0000) Subject: Move us a step further towards hp for continuous elements. X-Git-Tag: v8.0.0~11866 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=5e9bfa61b00ca97abc0c27251d334ce25fde3e09;p=dealii.git Move us a step further towards hp for continuous elements. git-svn-id: https://svn.dealii.org/trunk@12906 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/deal.II/source/dofs/hp_dof_handler.cc b/deal.II/deal.II/source/dofs/hp_dof_handler.cc index e052e7e56d..f810ceceff 100644 --- a/deal.II/deal.II/source/dofs/hp_dof_handler.cc +++ b/deal.II/deal.II/source/dofs/hp_dof_handler.cc @@ -2151,15 +2151,13 @@ namespace hp template <> void DoFHandler<2>::reserve_space () { + const unsigned int dim = 2; + Assert (finite_elements != 0, ExcNoFESelected()); Assert (finite_elements->size() > 0, ExcNoFESelected()); Assert (tria->n_levels() > 0, ExcInvalidTriangulation()); Assert (tria->n_levels() == levels.size (), ExcInternalError ()); - Assert (finite_elements->max_dofs_per_face() == 0, - ExcMessage ("hp finite elements are presently only supported " - "for discontinuous elements")); - // Release all space except the // active_fe_indices field which // we have to backup before @@ -2177,12 +2175,16 @@ namespace hp for (unsigned int level=0; leveln_levels(); ++level) { levels.push_back (new internal::hp::DoFLevel<2>); - std::swap (active_fe_backup[level], levels[level]->active_fe_indices); + std::swap (active_fe_backup[level], + levels[level]->active_fe_indices); } } + // QUAD (CELL) DOFs + // // count how much space we need - // on each level and set the + // on each level for the quad + // dofs and set the // dof_*_index_offset // data. initially set the latter // to an invalid index, and only @@ -2211,6 +2213,263 @@ namespace hp } + // LINE DOFS + // + // same here: count line dofs, + // then allocate as much space as + // we need and prime the linked + // list for lines (see the + // description in hp::DoFLevels) + // with the indices we will + // need. note that our task is + // more complicated since two + // adjacent cells may have + // different active_fe_indices, + // in which case we need to + // allocate *two* sets of line + // dofs for the same line + // + // the way we do things is that + // we loop over all active cells + // (these are the ones that have + // DoFs only anyway) and all + // their faces. We note in the + // user flags whether we have + // previously visited a face and + // if so skip it (consequently, + // we have to save and later + // restore the line flags) + { + std::vector saved_line_user_flags; + const_cast&>(*tria) + .save_user_flags_line (saved_line_user_flags); + const_cast&>(*tria).clear_user_flags_line (); + + // an array to hold how many + // slots (see the hp::DoFLevel + // class) we will have to store + // on each level + std::vector n_line_slots (tria->n_levels(), 0); + + for (active_cell_iterator cell=begin_active(); cell!=end(); ++cell) + for (unsigned int face=0; face::faces_per_cell; ++face) + if (! cell->face(face)->user_flag_set()) + { + // ok, face has not been + // visited. so we need to + // allocate space for + // it. let's see how much + // we need: we need one + // set if a) there is no + // neighbor behind this + // face, or b) the + // neighbor is not on the + // same level or further + // refined, or c) the + // neighbor is on the + // same level, but + // happens to have the + // same active_fe_index: + if (cell->at_boundary(face) + || + (cell->neighbor(face)->level() < cell->level()) + || + cell->neighbor(face)->has_children() + || + ((cell->neighbor(face)->level() == cell->level()) + && + !cell->neighbor(face)->has_children() + && + (cell->active_fe_index() == cell->neighbor(face)->active_fe_index()))) + // ok, one set of + // dofs. that makes + // one index, 1 times + // dofs_per_line + // dofs, and one stop + // index + n_line_slots[cell->level()] + += (*finite_elements)[cell->active_fe_index()].dofs_per_line + 2; + + // otherwise we do + // indeed need two + // sets, i.e. two + // indices, two sets of + // dofs, and one stop + // index: + else + n_line_slots[cell->level()] + += ((*finite_elements)[cell->active_fe_index()].dofs_per_line + + + (*finite_elements)[cell->neighbor(face)->active_fe_index()] + .dofs_per_line + + + 3); + + // mark this face as + // visited + cell->face(face)->set_user_flag (); + } + + // now that we know how many + // line dofs we will have to + // have on each level, allocate + // the memory. note that we + // allocate offsets for all + // lines, though only the + // active ones will have a + // non-invalid value later on + for (unsigned int level=0; leveln_levels(); ++level) + { + levels[level]->dof_line_index_offset + = std::vector (tria->n_raw_lines(level), + invalid_dof_index); + levels[level]->line_dofs + = std::vector (n_line_slots[level], + invalid_dof_index); + } + + // with the memory now + // allocated, loop over the + // cells again and prime the + // _offset values as well as + // the fe_index fields + const_cast&>(*tria).clear_user_flags_line (); + + std::vector next_free_line_slot (tria->n_levels(), 0); + + for (active_cell_iterator cell=begin_active(); cell!=end(); ++cell) + for (unsigned int face=0; face::faces_per_cell; ++face) + if (! cell->face(face)->user_flag_set()) + { + // same decision tree + // as before + if (cell->at_boundary(face) + || + (cell->neighbor(face)->level() < cell->level()) + || + cell->neighbor(face)->has_children() + || + ((cell->neighbor(face)->level() == cell->level()) + && + !cell->neighbor(face)->has_children() + && + (cell->active_fe_index() == cell->neighbor(face)->active_fe_index()))) + { + levels[cell->level()] + ->dof_line_index_offset[cell->face(face)->index()] + = next_free_line_slot[cell->level()]; + + // set first slot + // for this line to + // active_fe_index + // of this face + levels[cell->level()] + ->line_dofs[next_free_line_slot[cell->level()]] + = cell->active_fe_index(); + + // the next + // dofs_per_line + // indices remain + // unset for the + // moment (i.e. at + // invalid_dof_index). + // following this + // comes the stop + // index, which + // also is + // invalid_dof_index + // and therefore + // does not have to + // be explicitly + // set + + // finally, mark + // those slots as + // used + next_free_line_slot[cell->level()] + += (*finite_elements)[cell->active_fe_index()].dofs_per_line + 2; + } + else + { + levels[cell->level()] + ->dof_line_index_offset[cell->face(face)->index()] + = next_free_line_slot[cell->level()]; + + // set first slot + // for this line to + // active_fe_index + // of this face + levels[cell->level()] + ->line_dofs[next_free_line_slot[cell->level()]] + = cell->active_fe_index(); + + // the next + // dofs_per_line + // indices remain + // unset for the + // moment (i.e. at + // invalid_dof_index). + // + // then comes the + // fe_index for the + // neighboring + // cell: + levels[cell->level()] + ->line_dofs[next_free_line_slot[cell->level()] + + + (*finite_elements)[cell->active_fe_index()].dofs_per_line + + + 1] + = cell->neighbor(face)->active_fe_index(); + // then again a set + // of dofs that we + // need not set + // right now + // + // following this + // comes the stop + // index, which + // also is + // invalid_dof_index + // and therefore + // does not have to + // be explicitly + // set + + // finally, mark + // those slots as + // used + next_free_line_slot[cell->level()] + += ((*finite_elements)[cell->active_fe_index()].dofs_per_line + + + (*finite_elements)[cell->neighbor(face)->active_fe_index()] + .dofs_per_line + + + 3); + } + + // mark this face as + // visited + cell->face(face)->set_user_flag (); + } + + // we should have moved the + // cursor for each level to the + // total number of dofs on that + // level. check that + for (unsigned int level=0; leveln_levels(); ++level) + Assert (next_free_line_slot[level] == n_line_slots[level], + ExcInternalError()); + + // at the end, restore the user + // flags for the lines + const_cast&>(*tria) + .load_user_flags_line (saved_line_user_flags); + } + + +//TODO: allocate space for vertex dofs + // safety check: make sure that // the number of DoFs we // allocated is actually correct