From: wolf Date: Thu, 25 Jul 2002 16:49:15 +0000 (+0000) Subject: Do part of the work that was still missing. X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=4c3c531867e4fce66ca69cf21009f633f3362eee;p=dealii-svn.git Do part of the work that was still missing. git-svn-id: https://svn.dealii.org/trunk@6288 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/deal.II/source/fe/fe_nedelec.cc b/deal.II/deal.II/source/fe/fe_nedelec.cc index dea80725fa..a4570a157e 100644 --- a/deal.II/deal.II/source/fe/fe_nedelec.cc +++ b/deal.II/deal.II/source/fe/fe_nedelec.cc @@ -202,7 +202,7 @@ FE_Nedelec::clone() const template double FE_Nedelec::shape_value_component (const unsigned int i, - const Point &p, + const Point &p, const unsigned int component) const { Assert (idofs_per_cell, ExcIndexRange(i,0,this->dofs_per_cell)); @@ -462,30 +462,39 @@ FE_Nedelec::shape_grad_grad_component (const unsigned int i, template void FE_Nedelec::initialize_unit_support_points () { -//TODO: fix for higher orders. correct now for lowest order, all dimensions -// is this correct? all DoFs on lines, none on faces or bubbles? - - // all degrees of freedom are on - // edges, and their order is the - // same as the edges themselves - this->unit_support_points.resize(GeometryInfo::lines_per_cell * degree); - unsigned int index = 0; - for (unsigned int line=0; line::lines_per_cell; ++line) + switch (degree) { - const unsigned int - vertex_index_0 = GeometryInfo::vertices_adjacent_to_line(line,0), - vertex_index_1 = GeometryInfo::vertices_adjacent_to_line(line,1); - - const Point - vertex_0 = GeometryInfo::unit_cell_vertex(vertex_index_0), - vertex_1 = GeometryInfo::unit_cell_vertex(vertex_index_1); - - // place dofs equispaced - // between the vertices of each - // line - for (unsigned int d=0; dunit_support_points[index] - = (vertex_0*(d+1) + vertex_1*(degree-d)) / (degree+1); + case 1: + { + // all degrees of freedom are + // on edges, and their order + // is the same as the edges + // themselves + this->unit_support_points.resize(GeometryInfo::lines_per_cell); + for (unsigned int line=0; line::lines_per_cell; ++line) + { + const unsigned int + vertex_index_0 = GeometryInfo::vertices_adjacent_to_line(line,0), + vertex_index_1 = GeometryInfo::vertices_adjacent_to_line(line,1); + + const Point + vertex_0 = GeometryInfo::unit_cell_vertex(vertex_index_0), + vertex_1 = GeometryInfo::unit_cell_vertex(vertex_index_1); + + // place dofs right + // between the vertices + // of each line + this->unit_support_points[line] = (vertex_0 + vertex_1) / 2; + + break; + }; + }; + + default: + // no higher order + // elements implemented + // right now + Assert (false, ExcNotImplemented()); }; }; @@ -504,29 +513,37 @@ void FE_Nedelec<1>::initialize_unit_face_support_points () template void FE_Nedelec::initialize_unit_face_support_points () { -//TODO: fix for higher orders. correct now for lowest order, all dimensions -// is this correct? all DoFs on lines, none on faces or bubbles? - // do this the same as above, but - // for one dimension less - this->unit_face_support_points.resize(GeometryInfo::lines_per_cell * degree); - unsigned int index = 0; - for (unsigned int line=0; line::lines_per_cell; ++line) + switch (degree) { - const unsigned int - vertex_index_0 = GeometryInfo::vertices_adjacent_to_line(line,0), - vertex_index_1 = GeometryInfo::vertices_adjacent_to_line(line,1); + case 1: + { + // do this the same as above, but + // for one dimension less + this->unit_face_support_points.resize(GeometryInfo::lines_per_cell); + for (unsigned int line=0; line::lines_per_cell; ++line) + { + const unsigned int + vertex_index_0 = GeometryInfo::vertices_adjacent_to_line(line,0), + vertex_index_1 = GeometryInfo::vertices_adjacent_to_line(line,1); - const Point - vertex_0 = GeometryInfo::unit_cell_vertex(vertex_index_0), - vertex_1 = GeometryInfo::unit_cell_vertex(vertex_index_1); - - // place dofs equispaced - // between the vertices of each - // line - for (unsigned int d=0; dunit_face_support_points[index] - = (vertex_0*(d+1) + vertex_1*(degree-d)) / (degree+1); - }; + const Point + vertex_0 = GeometryInfo::unit_cell_vertex(vertex_index_0), + vertex_1 = GeometryInfo::unit_cell_vertex(vertex_index_1); + + // place dofs right + // between the vertices of each + // line + this->unit_face_support_points[line] = (vertex_0 + vertex_1) / 2; + }; + break; + }; + + default: + // no higher order + // elements implemented + // right now + Assert (false, ExcNotImplemented()); + }; }; @@ -535,15 +552,15 @@ template std::vector FE_Nedelec::get_dpo_vector(const unsigned int degree) { -//TODO: fix for higher orders. correct now for lowest order, all dimensions - std::vector dpo(dim+1, 0U); -// can this be done in a dimension independent and degree independent way? -// if DoFs are located only on lines, the the following is the correct way + Assert (degree == 1, ExcNotImplemented()); - // put all degrees of freedom on - // the lines, and in particular - // @p{degree} DoFs per line: + // for degree==1, put all degrees + // of freedom on the lines, and in + // particular @p{degree} DoFs per + // line: + std::vector dpo(dim+1, 0U); dpo[1] = degree; + return dpo; } @@ -551,15 +568,12 @@ FE_Nedelec::get_dpo_vector(const unsigned int degree) template UpdateFlags -FE_Nedelec::update_once (const UpdateFlags flags) const +FE_Nedelec::update_once (const UpdateFlags) const { -//TODO: think about what this actually means here??? - // for this kind of elements, only - // the values can be precomputed - // once and for all. set this flag - // if the values are requested at - // all - return (update_default | (flags & update_values)); + // even the values have to be + // computed on the real cell, so + // nothing can be done in advance + return update_default; } @@ -568,12 +582,12 @@ template UpdateFlags FE_Nedelec::update_each (const UpdateFlags flags) const { -//TODO: think about what this actually means here??? - UpdateFlags out = update_default; + if (flags & update_values) + out |= update_values | update_covariant_transformation; if (flags & update_gradients) - out |= update_gradients | update_covariant_transformation; + out |= update_gradients | update_covariant_transformation; if (flags & update_second_derivatives) out |= update_second_derivatives | update_covariant_transformation; @@ -588,78 +602,63 @@ FE_Nedelec::update_each (const UpdateFlags flags) const template typename Mapping::InternalDataBase * -FE_Nedelec::get_data (const UpdateFlags /*update_flags*/, - const Mapping &/*mapping*/, - const Quadrature &/*quadrature*/) const +FE_Nedelec::get_data (const UpdateFlags update_flags, + const Mapping &mapping, + const Quadrature &quadrature) const { - return 0; -//TODO -// // generate a new data object and -// // initialize some fields -// InternalData* data = new InternalData; - -// // check what needs to be -// // initialized only once and what -// // on every cell/face/subface we -// // visit -// data->update_once = update_once(update_flags); -// data->update_each = update_each(update_flags); -// data->update_flags = data->update_once | data->update_each; - -// const UpdateFlags flags(data->update_flags); -// const unsigned int n_q_points = quadrature.n_quadrature_points; - -// // some scratch arrays -// std::vector values(0); -// std::vector > grads(0); -// std::vector > grad_grads(0); - -// // initialize fields only if really -// // necessary. otherwise, don't -// // allocate memory -// if (flags & update_values) -// { -// values.resize (dofs_per_cell); -// data->shape_values.resize(dofs_per_cell, -// std::vector(n_q_points)); -// } - -// if (flags & update_gradients) -// { -// grads.resize (dofs_per_cell); -// data->shape_gradients.resize(dofs_per_cell, -// std::vector >(n_q_points)); -// } - -// // if second derivatives through -// // finite differencing is required, -// // then initialize some objects for -// // that -// if (flags & update_second_derivatives) -// data->initialize_2nd (this, mapping, quadrature); - -// // next already fill those fields -// // of which we have information by -// // now. note that the shape -// // gradients are only those on the -// // unit cell, and need to be -// // transformed when visiting an -// // actual cell -// if (flags & (update_values | update_gradients)) -// for (unsigned int i=0; ishape_values[renumber[k]][i] = values[k]; + // generate a new data object and + // initialize some fields + InternalData* data = new InternalData; + + // check what needs to be + // initialized only once and what + // on every cell/face/subface we + // visit + data->update_once = update_once(update_flags); + data->update_each = update_each(update_flags); + data->update_flags = data->update_once | data->update_each; + + const UpdateFlags flags(data->update_flags); + const unsigned int n_q_points = quadrature.n_quadrature_points; + + // initialize fields only if really + // necessary. otherwise, don't + // allocate memory + if (flags & update_values) + data->shape_values.reinit (dofs_per_cell, n_q_points); + + if (flags & update_gradients) + data->shape_gradients.reinit (dofs_per_cell, n_q_points); + + // if second derivatives through + // finite differencing is required, + // then initialize some objects for + // that + if (flags & update_second_derivatives) + data->initialize_2nd (this, mapping, quadrature); + + // next already fill those fields + // of which we have information by + // now. note that the shape values + // and gradients are only those on + // the unit cell, and need to be + // transformed when visiting an + // actual cell + for (unsigned int i=0; ishape_values[i][q][c] + = shape_value_component(i,quadrature.point(q),c); -// if (flags & update_gradients) -// for (unsigned int k=0; kshape_gradients[renumber[k]][i] = grads[k]; -// } -// return data; + if (flags & update_gradients) + for (unsigned int c=0; cshape_gradients[i][q][c] + = shape_grad_component(i,quadrature.point(q),c); + } + + return data; } @@ -671,45 +670,80 @@ FE_Nedelec::get_data (const UpdateFlags /*update_flags*/, template void -FE_Nedelec::fill_fe_values (const Mapping &/*mapping*/, - const typename DoFHandler::cell_iterator &/*cell*/, - const Quadrature &/*quadrature*/, - typename Mapping::InternalDataBase &/*mapping_data*/, - typename Mapping::InternalDataBase &/*fedata*/, - FEValuesData &/*data*/) const +FE_Nedelec::fill_fe_values (const Mapping &mapping, + const typename DoFHandler::cell_iterator &cell, + const Quadrature &quadrature, + typename Mapping::InternalDataBase &mapping_data, + typename Mapping::InternalDataBase &fedata, + FEValuesData &data) const { -//TODO!! -// // convert data object to internal -// // data for this class. fails with -// // an exception if that is not -// // possible -// InternalData &fe_data = dynamic_cast (fedata); + // convert data object to internal + // data for this class. fails with + // an exception if that is not + // possible + InternalData &fe_data = dynamic_cast (fedata); + + // get the flags indicating the + // fields that have to be filled + const UpdateFlags flags(fe_data.current_update_flags()); + + const unsigned int n_q_points = quadrature.n_quadrature_points; + + // fill shape function + // values. these are vector-valued, + // so we have to transform + // them. since the output format + // (in data.shape_values) is a + // sequence of doubles (one for + // each non-zero shape function + // value, and for each quadrature + // point, rather than a sequence of + // small vectors, we have to use a + // number of conversions + if (flags & update_values) + { + std::vector > shape_values (n_q_points); + for (unsigned int k=0; k