dolfin team mailing list archive

[Anders Logg <logg@xxxxxxxxx>]

Applied. Nice and clean demo.

-- 
Anders


On Thu, Apr 02, 2009 at 11:37:58AM -0600, Bartosz Sawicki wrote:
> Marie Rognes wrote:
>> Bartosz Sawicki wrote:
>>>>> I've written simple demo of curl-curl equation. It is based on  
>>>>> eddy-currents phenomena in low conducting materials. So far, in 
>>>>> the demo directory, there are no examples of this kind of PDE, 
>>>>> where Nedelec elements are used.
>>>>> If you found it useful, please consider placing it in the repository.
>>>>>
>>>>> The patch which I'm sending, contains cpp and python code. Cpp 
>>>>> works fine, but demo.py gives wrong results. I have no idea, 
>>>>> whats wrong in the python code. Please, take a look, this is 
>>>>> probably bug.
>>>>>
>>> >
>>>> (1) There seems to be a minus sign discrepancy between the two 
>>>> demos in the EddyCurrents part.
>>>
>>> You are right, there should be minus in both demos.
>>>
>>>
>>>> (2) In the python demo, replace
>>>>
>>>>    bc = DirichletBC(P1, zero, DirichletBoundary())
>>>>
>>>> by    bc = DirichletBC(PN, zero, DirichletBoundary())
>>>
>>> Oh, good! How could I miss that. This is the bug which I was looking  
>>> for. Now both examples works perfectly. Thank you very much.
>>>
>>> Do you want me to send improved patch? Or you can change this two  
>>> characters?
>>>
>>
>> I have no pushing privileges for dolfin -- so I guess someone higher up 
>> the hierarchy will get it done when they have time.
>> But, it is probably easiest  if you fix the bugs and send the perfect  
>> patch.
>
> The perfect patch is in the attachment.
>
> cheers,
>

> diff -r 124a218df971 demo/pde/curl-curl/cpp/CurrentDensity.form
> --- /dev/null	Thu Jan 01 00:00:00 1970 +0000
> +++ b/demo/pde/curl-curl/cpp/CurrentDensity.form	Thu Apr 02 11:34:05 2009 -0600
> @@ -0,0 +1,20 @@
> +# Copyright (C) 2009 Bartosz Sawicki.
> +# Licensed under the GNU LGPL Version 2.1.
> +#
> +# First added:  2009-03-30
> +# Last changed: 2009-03-30
> +#
> +# The bilinear form a(v, u) and linear form L(v) for
> +# curl operator.
> +#
> +# Compile this form with FFC: ffc -l dolfin CurrentDensity.form
> +
> +eL = VectorElement("CG", "tetrahedron", 1)
> +eN = FiniteElement("Nedelec", "tetrahedron", 0)
> +
> +v = TestFunction(eL)
> +u = TrialFunction(eL)
> +T = Function(eN)
> +
> +a = dot(v,u)*dx
> +L = dot(v,curl(T))*dx
> diff -r 124a218df971 demo/pde/curl-curl/cpp/EddyCurrents.form
> --- /dev/null	Thu Jan 01 00:00:00 1970 +0000
> +++ b/demo/pde/curl-curl/cpp/EddyCurrents.form	Thu Apr 02 11:34:05 2009 -0600
> @@ -0,0 +1,20 @@
> +# Copyright (C) 2009 Bartosz Sawicki.
> +# Licensed under the GNU LGPL Version 2.1.
> +#
> +# First added:  2009-03-30
> +# Last changed: 2009-03-30
> +#
> +# The bilinear form a(v, u) and linear form L(v) for
> +# curl-curl eddy currents equation.
> +#
> +# Compile this form with FFC: ffc -l dolfin EddyCurrents.form
> +
> +eN = FiniteElement("Nedelec", "tetrahedron", 0)
> +eL = VectorElement("Lagrange", "tetrahedron", 1)
> +
> +v = TestFunction(eN)
> +u = TrialFunction(eN)
> +dBdt = Function(eL)
> +
> +a = dot(rot(v), rot(u))*dx 
> +L = -dot(v, dBdt)*dx 
> diff -r 124a218df971 demo/pde/curl-curl/cpp/SConstruct
> --- /dev/null	Thu Jan 01 00:00:00 1970 +0000
> +++ b/demo/pde/curl-curl/cpp/SConstruct	Thu Apr 02 11:34:05 2009 -0600
> @@ -0,0 +1,16 @@
> +# This is an example of a simple SConstruct file for building
> +# programs against DOLFIN. To build this demo, just type 'scons'.
> +
> +import os, commands
> +
> +# Get compiler from pkg-config
> +compiler = commands.getoutput('pkg-config --variable=compiler dolfin')
> +
> +# Create a SCons Environment based on the main os environment
> +env = Environment(ENV=os.environ, CXX=compiler)
> +
> +# Get compiler flags from pkg-config
> +env.ParseConfig('pkg-config --cflags --libs dolfin')
> +
> +# Program name
> +env.Program('demo', 'main.cpp')
> diff -r 124a218df971 demo/pde/curl-curl/cpp/main.cpp
> --- /dev/null	Thu Jan 01 00:00:00 1970 +0000
> +++ b/demo/pde/curl-curl/cpp/main.cpp	Thu Apr 02 11:34:05 2009 -0600
> @@ -0,0 +1,98 @@
> +// Copyright (C) 2009 Bartosz Sawicki.
> +// Licensed under the GNU LGPL Version 2.1.
> +//
> +// First added:  2009-03-30
> +// Last changed: 2009-03-30
> +//
> +// Eddy currents phenomena in low conducting body can be 
> +// described using electric vector potential and curl-curl operator:
> +//    \nabla \times \nabla \times T = - \frac{\partial B}{\partial t}
> +// Electric vector potential defined as:
> +//    \nabla \times T = J
> +//
> +// Boundary condition 
> +//    J_n = 0, 
> +//    T_t=T_w=0, \frac{\partial T_n}{\partial n} = 0
> +// which is naturaly fulfilled for zero Dirichlet BC with Nedelec (edge) 
> +// elements.
> +
> +#include <dolfin.h>
> +#include "EddyCurrents.h"
> +#include "CurrentDensity.h"
> +  
> +using namespace dolfin;
> +
> +int main()
> +{
> +  
> +  // Homogenous external magnetic field (dB/dt)
> +  class Source : public Function
> +  {
> +  public:
> +    void eval(real* values, const double* x) const
> +    {
> +      values[0] = 0.0;
> +      values[1] = 0.0;
> +      values[2] = 1.0;
> +    }
> +  };
> +
> +  // Zero Dirichlet BC
> +  class Zero : public Function
> +  {
> +  public:
> +    void eval(real* values, const double* x) const
> +    {
> +      values[0] = 0.0;
> +      values[1] = 0.0;
> +      values[2] = 0.0;
> +    }
> +  };
> +
> +  // Everywhere on external surface
> +  class DirichletBoundary: public SubDomain
> +  {
> +    bool inside(const real* x, bool on_boundary) const
> +    {
> +      return on_boundary;
> +    }
> +  };
> +
> +  // Create demo mesh
> +  UnitSphere mesh(10);
> +
> +  // Define functions
> +  Source dBdt;
> +  Zero zero;
> +  Function T;
> +  Function J;
> +
> +  // Define function space and boundary condition
> +  EddyCurrentsFunctionSpace V(mesh);
> +  DirichletBoundary boundary;
> +  DirichletBC bc(V, zero, boundary);
> +
> +  // Use forms to define variational problem
> +  EddyCurrentsBilinearForm a(V,V);
> +  EddyCurrentsLinearForm L(V);
> +  L.dBdt = dBdt;
> +  VariationalProblem problem (a, L,  bc);
> +  
> +  // Solve problem using default solver
> +  problem.solve(T);
> +
> +  // Define variational problem for current density (J)
> +  CurrentDensityFunctionSpace V1(mesh);
> +  CurrentDensityBilinearForm a1(V1,V1);
> +  CurrentDensityLinearForm L1(V1);
> +  L1.T = T;
> +  VariationalProblem problem1(a1, L1);
> +
> +  // Solve problem using default solver
> +  problem1.solve(J);
> +
> +  // Plot solution
> +  plot(J);
> +
> +  return 0;
> +}
> diff -r 124a218df971 demo/pde/curl-curl/python/demo.py
> --- /dev/null	Thu Jan 01 00:00:00 1970 +0000
> +++ b/demo/pde/curl-curl/python/demo.py	Thu Apr 02 11:34:05 2009 -0600
> @@ -0,0 +1,64 @@
> +""" Eddy currents phenomena in low conducting body can be 
> +described using electric vector potential and curl-curl operator:
> +   \nabla \times \nabla \times T = - \frac{\partial B}{\partial t}
> +Electric vector potential defined as:
> +   \nabla \times T = J
> +
> +Boundary condition:
> +   J_n = 0, 
> +   T_t=T_w=0, \frac{\partial T_n}{\partial n} = 0
> +which is naturaly fulfilled for zero Dirichlet BC with Nedelec (edge) 
> +elements.
> +"""
> +
> +__author__ = "Bartosz Sawicki (sawickib@xxxxxxxxxxxxx)"
> +__date__ = "2009-04-02"
> +__copyright__ = "Copyright (C) 2009 Bartosz Sawicki"
> +__license__  = "GNU LGPL Version 2.1"
> +
> +from dolfin import *
> +
> +# Create mesh
> +mesh = UnitSphere(10)
> +
> +# Define function spaces
> +PN = FunctionSpace(mesh, "Nedelec", 0)
> +P1 = VectorFunctionSpace(mesh, "CG", 1)
> +
> +# Define test and trial functions
> +v0 = TestFunction(PN)
> +u0 = TrialFunction(PN)
> +v1 = TestFunction(P1)
> +u1 = TrialFunction(P1)
> +
> +# Define functions
> +dBdt = Function(P1, ("0.0", "0.0", "1.0"))
> +zero = Function(P1, ("0.0", "0.0", "0.0"))
> +T = Function(PN)
> +J = Function(P1)
> +
> +# Dirichlet boundary
> +class DirichletBoundary(SubDomain):
> +    def inside(self, x, on_boundary):
> +        return on_boundary
> +
> +# Boundary condition
> +bc = DirichletBC(PN, zero, DirichletBoundary())
> +
> +# Eddy currents equation (using potential T) 
> +Teqn = (dot(rot(v0), rot(u0))*dx, -dot(v0, dBdt)*dx)
> +Tproblem = VariationalProblem( Teqn[0], Teqn[1], bc)
> +T = Tproblem.solve()
> +
> +# Current density equation
> +Jeqn = (dot(v1, u1)*dx, dot(v1, curl(T))*dx)
> +Jproblem = VariationalProblem( Jeqn[0], Jeqn[1])
> +J = Jproblem.solve()
> +
> +# Plot solution
> +plot(J)
> +
> +# Hold plot
> +interactive()
> +
> +

> _______________________________________________
> DOLFIN-dev mailing list
> DOLFIN-dev@xxxxxxxxxx
> http://www.fenics.org/mailman/listinfo/dolfin-dev

Attachment: signature.asc
Description: Digital signature