dolfin team mailing list archive

[Anders Logg <logg@xxxxxxxxx>]

On Mon, Aug 11, 2008 at 06:21:05PM +0200, Dag Lindbo wrote:
> 
> 
> Anders Logg wrote:
> > On Mon, Jun 30, 2008 at 11:44:40AM +0200, Dag Lindbo wrote:
> >> Anders Logg wrote:
> >>> On Tue, Jun 17, 2008 at 10:52:06PM +0200, Dag Lindbo wrote:
> >>>> Hello again!
> >>>>
> >>>> Is it correct that Periodic BCs won't work for vector-valued systems?
> >>>> There is a comment to suggest this in PeriodicBC.cpp. Is there a
> >>>> reasonable workaround for the needy?
> >>>>
> >>>> /Dag
> >>> I don't know of an easy workaround.
> >>>
> >>> But maybe you can take a look at the code in PeriodicBC.cpp to see if
> >>> it's possible to understand what it does and then extend it? I guess
> >>> the extension will be nontrivial, but the code is fairly well
> >>> documented so it should not be impossible.
> >> Sorry for the basic question, but I'm struggling to figure out the role
> >> SubSystems play in PeriodicBCs.
> >>
> >> E.g. on a Taylor-Hood element, SubSystem(1) corresponds to pressure - a
> >> scalar. I expected that it would be possible to apply the PBC to this
> >> subsystem. But if I do
> >> PeriodicBC bc(mesh, sub_domain, SubSystem(1)),
> >> and try to apply it I still get an abort that there are "more than one
> >> dof associated with coordinate. Did you forget to specify the subsystem?".
> >>
> >> I would greatly appreciate if someone could help to clarify this.
> >>
> >> /Dag
> > 
> > Is this still an issue?
> > 
> 
> Yes, I have not been able to make more sense of this. Sorry. /Dag

I found the following comment in PeriodicBC:

  // FIXME: Make this work for non-scalar subsystems, like vector-valued
  // FIXME: Lagrange where more than one per element is associated with
  // FIXME: each coordinate. Note that globally there may very well be
  // FIXME: more than one dof per coordinate (for conforming elements).

This seems to be related to the error message that you get.

What is done in PeriodicBC is that each coordinate is mapped to a
degree of freedom. This way, we may find if the coordinate of one dof
is mapped to that same dof and set them equal. I think it should be
possible to extend to the case where there are two dofs or more
associated with each coordinate (like for vector elements). We just
need to include the component/subsystem in the key.

Anyway, I've managed to use PeriodicBC before for setting periodic
boundary conditions for a complex-valued Stokes problem. I've attached
part of the code so you can compare. Tell me if you need more.

-- 
Anders

#include <map>

#include <dolfin.h>
#include "ComplexStokes.h"
#include "Average.h"

using namespace dolfin;

// Antal obekanta:
//
// Hastighet: 2*2*(v + e) 
// Tryck:     2*v
// Totalt:    6v + 4e
//
// v = (n+1)*(n+1)
// e = 2*(n+1)*n + n*n = n*(3n + 2)
//
// Totalt: 6*(n+1)*(n+1) + 4*n*(3*n + 2)
//
//
//     10    2006
//     20    7606
//     30   16806
//     40   29606
//     50   46006
//     60   66006
//     70   89606
//     80  116806
//     90  147606
//    100  182006
//    110  220006
//    120  261606
//    130  306806
//    140  355606
//    150  408006

class NoslipBoundary : public SubDomain
{
public:
  
  NoslipBoundary(real b) : b(b) {}

  bool inside(const real* x, bool on_boundary) const
  {
    return ( (x[1] < DOLFIN_EPS && x[0] < b) ||
             (x[1] < DOLFIN_EPS && x[0] > (1.0 - b)) ||
             (x[1] > (1.0 - DOLFIN_EPS) && x[0] < b) ||
             (x[1] > (1.0 - DOLFIN_EPS) && x[0] > (1.0 - b)) ||
             (x[0] < DOLFIN_EPS && x[1] < b) ||
             (x[0] < DOLFIN_EPS && x[1] > (1.0 - b)) ||
             (x[0] > (1.0 - DOLFIN_EPS) && x[1] < b) ||
             (x[0] > (1.0 - DOLFIN_EPS) && x[1] > (1.0 - b)) ) && on_boundary;
  }

private:

  real b;
  
};

class NoslipBoundaryY : public SubDomain
{
public:

  NoslipBoundaryY(real b) : b(b) {}

  bool inside(const real* x, bool on_boundary) const
  {
    return ( (x[1] < DOLFIN_EPS && x[0] > b && x[0] < (1.0 - b)) ||
             (x[1] > (1.0 - DOLFIN_EPS) && x[0] > b && x[0] < (1.0 - b)) );
  }

private:

  real b;
  
};

class PeriodicBoundary : public SubDomain
{
public:

  PeriodicBoundary(real b) : b(b) {}
  
  bool inside(const real* x, bool on_boundary) const
  {
    return ( (x[0] < DOLFIN_EPS && x[1] > b && x[1] < (1.0 - b)) ||
             (x[0] > b && x[0] < (1.0 - b) && x[1] < DOLFIN_EPS) ) && on_boundary;
  }

  void map(const real* x, real* y) const
  {
    if (x[0] > (1.0 - DOLFIN_EPS))
    {
      // Map from right to left
      y[0] = x[0] - 1.0;
      y[1] = x[1];
    }
    else
    {
      // Map from up to down
      y[0] = x[0];
      y[1] = x[1] - 1.0;
    }
  }

private:

  real b;

};

class PressureBoundary : public SubDomain
{
public:
  
  PressureBoundary(real h) : h(h) {}
  
  bool inside(const real* x, bool on_boundary) const
  {
    real w = 2.0*h;
    
    return x[0] < w && x[1] < w;
  }
  
private:
  
  real h;
  
};

class NoslipFunction : public Function
{
public:
  
  NoslipFunction(Mesh& mesh) : Function(mesh) {}
  
  void eval(real* values, const real* x) const
  {
    values[0] = 0.0;
    values[1] = 0.0;
  }
  
};

class NoslipFunctionY : public Function
{
public:
  
  NoslipFunctionY(Mesh& mesh) : Function(mesh) {}
  
  void eval(real* values, const real* x) const
  {
    values[0] = 0.0;
  }
  
};

class ForceFunctionReal : public Function
{
public:
  
  ForceFunctionReal(Mesh& mesh) : Function(mesh) {}
  
  void eval(real* values, const real* x) const
  {
    values[0] = 1.0;
    values[1] = 0.0;
  }
  
};

class ForceFunctionImag : public Function
{
public:
  
  ForceFunctionImag(Mesh& mesh) : Function(mesh) {}
  
  void eval(real* values, const real* x) const
  {
    values[0] = 0.0;
    values[1] = 0.0;
  }
  
};

void solve(real aa, real mu, real rho, real omega)
{
  // Width of wall
  real bb = 0.5 - aa;

  // Create mesh
  unsigned int n = 128;
  real h = 1.0 / static_cast<real>(n);
  UnitSquare mesh(n, n);

  // Create sub domains for boundary conditions
  NoslipBoundary noslip_boundary(bb);
  NoslipBoundaryY noslip_boundary_y(bb);
  PressureBoundary pressure_boundary(h);
  PeriodicBoundary periodic_boundary(bb);

  // Create functions for boundary conditions
  NoslipFunction noslip_function(mesh);
  NoslipFunctionY noslip_function_y(mesh);
  Function zero(mesh, 0.0);

  // Define sub systems for boundary conditions
  SubSystem velocity_r(0);
  SubSystem velocity_i(1);
  SubSystem pressure_r(2);
  SubSystem pressure_i(3);
  SubSystem velocity_rx(0, 0);
  SubSystem velocity_ry(0, 1);
  SubSystem velocity_ix(1, 0);
  SubSystem velocity_iy(1, 1);

  // No-slip boundary condition for velocity
  DirichletBC bc0(noslip_function, mesh, noslip_boundary, velocity_r);
  DirichletBC bc1(noslip_function, mesh, noslip_boundary, velocity_i);
  DirichletBC bc10(noslip_function_y, mesh, noslip_boundary_y, velocity_ry);
  DirichletBC bc11(noslip_function_y, mesh, noslip_boundary_y, velocity_iy);

  // Pressure boundary condition
  DirichletBC bc2(zero, mesh, pressure_boundary, pressure_r);
  DirichletBC bc3(zero, mesh, pressure_boundary, pressure_i);

  // Periodic boundary conditions
  PeriodicBC bc4(mesh, periodic_boundary, velocity_rx);
  PeriodicBC bc5(mesh, periodic_boundary, velocity_ry);
  PeriodicBC bc6(mesh, periodic_boundary, velocity_ix);
  PeriodicBC bc7(mesh, periodic_boundary, velocity_iy);
  PeriodicBC bc8(mesh, periodic_boundary, pressure_r);
  PeriodicBC bc9(mesh, periodic_boundary, pressure_i);

  // Create force terms
  ForceFunctionReal fr(mesh);
  ForceFunctionImag fi(mesh);
  
  // Create constants
  Function mu_function(mesh, mu);
  Function rho_function(mesh, rho);
  
  // Set up variational problem
  ComplexStokesBilinearForm a(mu_function, rho_function);
  ComplexStokesLinearForm L(fr, fi);

  // Assemble linear system
  Matrix A;
  Vector b;
  assemble(A, a, mesh);
  assemble(b, L, mesh);
  
  // Apply boundary conditions
  bc0.apply(A, b, a);
  bc1.apply(A, b, a);
  //bc2.apply(A, b, a);
  //bc3.apply(A, b, a);
  bc4.apply(A, b, a);
  bc5.apply(A, b, a);
  bc6.apply(A, b, a);
  bc7.apply(A, b, a);
  bc8.apply(A, b, a);
  bc9.apply(A, b, a);
  bc10.apply(A, b, a);
  bc11.apply(A, b, a);
  
  // Solve linear system
  Vector x;
  solve(A, x, b);
  Function u(mesh, x, a);

  // Extract sub functions
  Function wr = u[0];
  Function wi = u[1];
  Function pr = u[2];
  Function pi = u[3];

  // Compute average of x-components of wr and wi
  AverageFunctional Mr(wr);
  AverageFunctional Mi(wi);
  real Wr = assemble(Mr, mesh);
  real Wi = assemble(Mi, mesh);
  message("Result: %.15e %.15e %.15e %.15e", aa, omega, Wr, Wi);

  /*

  // Plot solution
  plot(wr);
  plot(wi);
  plot(pr);
  plot(pi);
  
  // Save matrix
  //File file("matrix.m");
  //file << A;
  
  // Save solution to XML
  File f0("velocity_real.xml");
  f0 << wr;

  File f1("velocity_imag.xml");
  f1 << wi;

  File f2("pressure_real.xml");
  f2 << pr;

  File f3("pressure_imag.xml");
  f3 << pi;

  // Save solution to PVD
  File f4("velocity_real.pvd");
  f4 << wr;

  File f5("velocity_imag.pvd");
  f5 << wi;

  File f6("pressure_real.pvd");
  f6 << pr;

  File f7("pressure_imag.pvd");
  f7 << pi;

  */
}

int main(int argc, char* argv[])
{
  if (argc != 3)
    error("Usage: stokes a omega");

  real a = atof(argv[1]);
  real omega = atof(argv[2]);
  message("Solving for a = %g, omega = %g", a, omega);
            
  // Parameters
  real rho_s = 1.0e3;
  real rho_f = 1.3;
  real mu    = 1.7e-5;
  real h     = 1.0e-4;

  // Scaled parameters
  real mu_tilde = mu / (h*h*omega*rho_f);
  real rho_hat  = rho_f / rho_s;
  
  solve(a, mu_tilde, rho_hat, omega);
}

Attachment: signature.asc
Description: Digital signature