Dear John, I’m not sure that there is an interpolation function that would work in that way for Nedelec elements (there is VectorTools::project_boundary_values_curl_conforming_l2() <https://dealii.org/current/doxygen/deal.II/group__constraints.html#gacdba6b0a8a06f9c9ceff9f11ef05a203> for Neumann boundary conditions, not that this is particularly useful for you here). As always we’d be happy to accept a contribution to the library for the extension that you require!
You don’t mention specifically what you’re trying to model, or where this interpolated function is to be used, so I’m just going to throw out a couple of things for you to think about and decide if they’re valid or of use to your specific application. What I’ve done in the past is to use a dim-component Function for the free current source, and I’d written a check (which I’ve pasted below) to verify that I’d not set up this source function incorrectly. Admittedly, I had a relatively straight-forward geometry so it wasn’t challenging to write such a source function. If you’ve got a more complex setup with the current source derived from some secondary finite element problem then maybe such an approach is not appropriate. But if it so happens that you can define or solve for a scalar potential whose gradient is the source current, then this is divergence free if the scalar potential is approximated using Lagrange polynomials (i.e. FE_Q <https://dealii.org/current/doxygen/deal.II/classFE__Q.html> . You can then use an FEFieldFunction <https://dealii.org/current/doxygen/deal.II/classFunctions_1_1FEFieldFunction.html> to compute this gradient at each point within the current source, to act as a source term for the magnetic vector potential problem. Best, Jean-Paul template <int dim> void WireMagneticField<dim>::verify_source_terms () const { TimerOutput::Scope timer_scope (computing_timer, "Verify source terms"); hp::FEFaceValues<dim> hp_fe_face_values (mapping_collection, fe_collection, qf_collection_face, update_quadrature_points | update_normal_vectors | update_JxW_values); // To verify that the source terms are divergence free, // we make use of the divergence theorem: // 0 = \int_{\Omega_{e}} div (J_f) // = \int_{\partial\Omega_{e}} J_f . n typename hp::DoFHandler<dim>::active_cell_iterator cell = hp_dof_handler.begin_active(), endc = hp_dof_handler.end(); for (; cell!=endc; ++cell) { if (cell->subdomain_id() != this_mpi_process) continue; double div_J_f = 0.0; for (unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face) { hp_fe_face_values.reinit(cell,face); const FEFaceValues<dim> &fe_face_values = hp_fe_face_values.get_present_fe_values(); const unsigned int &n_fq_points = fe_face_values.n_quadrature_points; std::vector<Vector<double> > source_values (n_fq_points, Vector<double>(dim)); function_free_current_density.vector_value_list (fe_face_values.get_quadrature_points(), source_values); for (unsigned int fq_point=0; fq_point<n_fq_points; ++fq_point) { const Tensor<1,dim> J_f ({source_values[fq_point][0], source_values[fq_point][1], source_values[fq_point][2]}); // Note: J_f must be divergence free! const double JxW = fe_face_values.JxW(fq_point); const Tensor<1,dim> &N = fe_face_values.normal_vector(fq_point); div_J_f += (J_f*N) * JxW; } } AssertThrow(std::abs(div_J_f) < 1e-9, ExcMessage("Source term is not divergence free!")); } } > On 16. Apr 2021, at 19:36, John Smith <[email protected]> wrote: > > Dear Simon, > > Thank you for your reply. Your suggestion definitely works. I used functions > that output one component when I was interpolating scalar potentials. It > works well. > > However, the FE_Nedelec is different. It implements edge elements. They are > vector-based. That is, functions are represented by a superposition of > vector-valued shape functions: > \vec{A} = \sum u_i vec{N}_i . > Therefore, the output of the "value" method in “class CustomFunction : public > Function” must be vector-valued. Three components in a three-dimensional > space. Otherwise, there is no point in interpolation. > > This kind of vectorial approximations is a bread-and-butter topic in > magnetics. See, for example, equation (7) in: > > https://ieeexplore.ieee.org/document/497322 > <https://ieeexplore.ieee.org/document/497322> > > In short, the source, i.e the current vector potential, must be projected on > the space spanned by the vector-valued shape functions. Otherwise, the > simulation is numerically unstable. The last, from my experience, is > definitely true. > > Best, > John > > On Fri, Apr 16, 2021 at 5:28 PM [email protected] > <mailto:[email protected]> <[email protected] > <mailto:[email protected]>> wrote: > Hi, > I don't have any experience with FE_Nedelec, but the error message states > that the function you are trying to interpolate has 1 component, while the > element has dim components. > > If you implemented your own Function, you need to make sure it has > dim-components by calling the constructor of Function with dim as inargument: > > template <int dim> > class CustomFunction : public Function<dim> > { > public: > CustomFunction() > : Function<dim>(dim) // The function has dim components > {} > > double > value(const Point<dim> & point, > const unsigned int component = 0) const override > { > // Return value will depend on component here > return 0; > } > }; > > > Best, > Simon > > On Friday, April 16, 2021 at 5:07:01 PM UTC+2 [email protected] > <mailto:[email protected]> wrote: > Hello, > > > It seems I am unable to find a function similar to VectorTools::interpolate > for FE_Nedelec finite elements. The existing implementation of this functions > gives the following run-time error: > > > > An error occurred in line <556> of file > </home/john/dealii-9.1.1/include/deal.II > > /numerics/vector_tools.templates.h> in function > > void dealii::VectorTools::interpolate(const dealii::Mapping<dim, spacedim>&, > > const DoFHandlerType<dim, spacedim>&, const dealii::Function<spacedim, > typename > > VectorType::value_type>&, VectorType&, const dealii::ComponentMask&) [with > int > > dim = 3; int spacedim = 3; VectorType = dealii::Vector<double>; > DoFHandlerType = > > dealii::DoFHandler; typename VectorType::value_type = double] > > The violated condition was: > > dof_handler.get_fe().n_components() == function.n_components > > Additional information: > > Dimension 3 not equal to 1. > > > > > Looks like VectorTools::interpolate does not like the vector-valued finite > elements. Could you please point me in the right direction? > > > > Best, > > John > > > -- > The deal.II project is located at http://www.dealii.org/ > <http://www.dealii.org/> > For mailing list/forum options, see > https://groups.google.com/d/forum/dealii?hl=en > <https://groups.google.com/d/forum/dealii?hl=en> > --- > You received this message because you are subscribed to a topic in the Google > Groups "deal.II User Group" group. > To unsubscribe from this topic, visit > https://groups.google.com/d/topic/dealii/wG3j4yD1rpw/unsubscribe > <https://groups.google.com/d/topic/dealii/wG3j4yD1rpw/unsubscribe>. > To unsubscribe from this group and all its topics, send an email to > [email protected] > <mailto:[email protected]>. > To view this discussion on the web visit > https://groups.google.com/d/msgid/dealii/6d051c9e-8840-4e5d-8021-f8007a1dc543n%40googlegroups.com > > <https://groups.google.com/d/msgid/dealii/6d051c9e-8840-4e5d-8021-f8007a1dc543n%40googlegroups.com?utm_medium=email&utm_source=footer>. > > -- > The deal.II project is located at http://www.dealii.org/ > <http://www.dealii.org/> > For mailing list/forum options, see > https://groups.google.com/d/forum/dealii?hl=en > <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] > <mailto:[email protected]>. > To view this discussion on the web visit > https://groups.google.com/d/msgid/dealii/CAH3L9AVtN%3Dtaqi2YiTjfoyfhxo-%2B9XnNzpBqNXSHqk5wX4iUtQ%40mail.gmail.com > > <https://groups.google.com/d/msgid/dealii/CAH3L9AVtN%3Dtaqi2YiTjfoyfhxo-%2B9XnNzpBqNXSHqk5wX4iUtQ%40mail.gmail.com?utm_medium=email&utm_source=footer>. -- 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/A3354A92-DBBF-47ED-A75A-239AA49598E5%40gmail.com.
