I was reading the documentation and came across CompiledSubdomain here:
https://fenicsproject.org/pub/tutorial/html/._ftut1013.html#___sec96
Is there a way I can implement Periodic Boundaries on a simple rectangular mesh using CompiledSubdomain ? I looked for it, but did not find anything.
Hi,
You can implement periodic boundary conditions using the SubDomain class (specifically its map method) as described in Periodic homogenization of linear elastic materials — Numerical tours of continuum mechanics using FEniCS master documentation
@bhaveshshrimali Yes, for now I have implemented PBC using SubDomain. The link
I have given in my original post mentions that CompiledSubDomain eliminates the need of python calls from C++, which is why I wanted to used it. Given examples however cover only simple Dirichlet expressions.
I see. If performance is a significant bottleneck you can directly use the underlying C++ class. See CompiledSubdomain using C++ class in case if you haven’t already. You can view the complete list of methods at Bitbucket. Though map and inside are the only relevant methods that you may need.
Edit: I haven’t done that myself yet. But should be straightforward given that you have already gotten the idea correct. If I get some time, I can make a simple MWE.
@bhaveshshrimali So I managed to cook up the following by looking at committor/periodic_boundary_conditions.py at 8b05d258c538e57385f4be95c196be2d2390f4bd · jmbr/committor · GitHub
import fenics as fn
CPP_CODE = '''
class PeriodicBoundary : public SubDomain {
public:
bool inside(const Array<double>& x, bool on_boundary) const {
bool left = std::fabs(x[0]) < DOLFIN_EPS;
bool bottom = std::fabs(x[1]) < DOLFIN_EPS;
bool right = std::fabs(x[0]-${lengthx}) < DOLFIN_EPS;
bool top = std::fabs(x[1]-${lengthy}) < DOLFIN_EPS;
return ( left || bottom ) && (!( right || top )) && on_boundary;
}
void map(const Array<double>& x, Array<double>& y) const {
y[0] = x[0];
y[1] = x[1];
if (std::fabs(x[0]- (${lengthx})) < DOLFIN_EPS)
y[0] -= ${lengthx};
if (std::fabs(x[1]- (${lengthy})) < DOLFIN_EPS)
y[1] -= ${lengthy};
}
};
'''
# Compiled periodic boundary condition
def compiled_pbc(lx,ly):
import string
template = string.Template(CPP_CODE)
code = template.safe_substitute(lengthx=f'{lx:.5E}',lengthy=f'{lx:.5E}')
print(code)
return fn.CompiledSubDomain(code)
But I end up with the following C++ error. Can you help me with this? I only know C.
------------------- Start compiler output ------------------------
/tmp/tmpql5fnaia/dolfin_subdomain_c6bd6ce80bdfde059f2a510d46ad0e92.cpp: In member function ‘virtual bool dolfin::dolfin_subdomain_c6bd6ce80bdfde059f2a510d46ad0e92::inside(Eigen::Ref<const Eigen::Matrix<double, -1, 1> >, bool) const’:
/tmp/tmpql5fnaia/dolfin_subdomain_c6bd6ce80bdfde059f2a510d46ad0e92.cpp:64:1: error: expected primary-expression before ‘class’
class PeriodicBoundary : public SubDomain {
^~~~~
/tmp/tmpql5fnaia/dolfin_subdomain_c6bd6ce80bdfde059f2a510d46ad0e92.cpp:64:1: error: expected ‘;’ before ‘class’
------------------- End compiler output ------------------------Here’s the C++ snippet it tries to compile. Note how the method bool inside(const Eigen:...) is incomplete. Looking at the source code here, python class CompiledSubDomain only appears to support simple expressions for bool inside cpp function.
// Based on https://gcc.gnu.org/wiki/Visibility
#if defined _WIN32 || defined __CYGWIN__
#ifdef __GNUC__
#define DLL_EXPORT __attribute__ ((dllexport))
#else
#define DLL_EXPORT __declspec(dllexport)
#endif
#else
#define DLL_EXPORT __attribute__ ((visibility ("default")))
#endif
#include <dolfin/common/Array.h>
#include <dolfin/math/basic.h>
#include <dolfin/mesh/SubDomain.h>
#include <Eigen/Dense>
// cmath functions
using std::cos;
using std::sin;
using std::tan;
using std::acos;
using std::asin;
using std::atan;
using std::atan2;
using std::cosh;
using std::sinh;
using std::tanh;
using std::exp;
using std::frexp;
using std::ldexp;
using std::log;
using std::log10;
using std::modf;
using std::pow;
using std::sqrt;
using std::ceil;
using std::fabs;
using std::floor;
using std::fmod;
using std::max;
using std::min;
const double pi = DOLFIN_PI;
namespace dolfin
{
class dolfin_subdomain_c6bd6ce80bdfde059f2a510d46ad0e92 : public SubDomain
{
public:
dolfin_subdomain_c6bd6ce80bdfde059f2a510d46ad0e92()
{
}
// Return true for points inside the sub domain
bool inside(const Eigen::Ref<const Eigen::VectorXd> x, bool on_boundary) const final
{
return
class PeriodicBoundary : public SubDomain {
public:
bool inside(const Array<double>& x, bool on_boundary) const {
bool left = std::fabs(x[0]) < DOLFIN_EPS;
bool bottom = std::fabs(x[1]) < DOLFIN_EPS;
bool right = std::fabs(x[0]-1.00000E+00) < DOLFIN_EPS;
bool top = std::fabs(x[1]-1.00000E+00) < DOLFIN_EPS;
return ( left || bottom ) && (!( right || top )) && on_boundary;
}
void map(const Array<double>& x, Array<double>& y) const {
y[0] = x[0];
y[1] = x[1];
if (std::fabs(x[0]- (1.00000E+00)) < DOLFIN_EPS)
y[0] -= 1.00000E+00;
if (std::fabs(x[1]- (1.00000E+00)) < DOLFIN_EPS)
y[1] -= 1.00000E+00;
}
};
;
}
void set_property(std::string name, double value)
{
}
double get_property(std::string name) const
{
return 0.0;
}
};
}
extern "C" DLL_EXPORT dolfin::SubDomain * create_dolfin_subdomain_c6bd6ce80bdfde059f2a510d46ad0e92()
{
return new dolfin::dolfin_subdomain_c6bd6ce80bdfde059f2a510d46ad0e92;
}Since the CompiledSubdomain class only supports simple inside expressions, I constructed up a new CompiledPBC class to get periodic boundary conditions working with C++ for two dimensional rectangular domains of size [lx \times ly]. Here’s the source code. Put it in compiledpbc.py file to use.
@dokken @bhaveshshrimali Any tips on maybe simplifying it?
Usage:
from compiledpbc import CompiledPBC
pbc = CompiledPBC(lx,ly)
Source:
# -*- coding: utf-8 -*-
# Copyright (C) 2017 Chris N. Richardson and Garth N. Wells
#
# This file is part of DOLFIN.
#
# DOLFIN is free software: you can redistribute it and/or modify
# it under the terms of the GNU Lesser General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# DOLFIN is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU Lesser General Public License for more details.
#
# You should have received a copy of the GNU Lesser General Public License
# along with DOLFIN. If not, see <http://www.gnu.org/licenses/>.
import dolfin.cpp as cpp
from dolfin.cpp import MPI
from dolfin.cpp.log import log, LogLevel
from dolfin.jit.jit import compile_class, _math_header
def jit_generate(class_data, module_name, signature, parameters):
template_code = """
// Based on https://gcc.gnu.org/wiki/Visibility
#if defined _WIN32 || defined __CYGWIN__
#ifdef __GNUC__
#define DLL_EXPORT __attribute__ ((dllexport))
#else
#define DLL_EXPORT __declspec(dllexport)
#endif
#else
#define DLL_EXPORT __attribute__ ((visibility ("default")))
#endif
#include <dolfin/common/Array.h>
#include <dolfin/math/basic.h>
#include <dolfin/mesh/SubDomain.h>
#include <Eigen/Dense>
{math_header}
namespace dolfin
{{
class {classname} : public SubDomain
{{
public:
{members}
{classname}()
{{
{constructor}
}}
// Return true for points inside the sub domain
bool inside(const Eigen::Ref<const Eigen::VectorXd> x, bool on_boundary) const final
{{
bool left = std::fabs(x[0]) < DOLFIN_EPS;
bool bottom = std::fabs(x[1]) < DOLFIN_EPS;
bool right = std::fabs(x[0] - {lengthx}) < DOLFIN_EPS;
bool top = std::fabs(x[1] - {lengthy}) < DOLFIN_EPS;
return ( left || bottom ) && (!( right || top )) && on_boundary;
}}
void map(const Eigen::Ref<const Eigen::VectorXd> x, Eigen::Ref<Eigen::VectorXd> y) const
{{
y[0] = x[0];
y[1] = x[1];
if (std::fabs(x[0] - ({lengthx})) < DOLFIN_EPS)
y[0] -= {lengthx};
if (std::fabs(x[1] - ({lengthy})) < DOLFIN_EPS)
y[1] -= {lengthy};
}}
void set_property(std::string name, double value)
{{
{set_props}
}}
double get_property(std::string name) const
{{
{get_props}
return 0.0;
}}
}};
}}
extern "C" DLL_EXPORT dolfin::SubDomain * create_{classname}()
{{
return new dolfin::{classname};
}}
"""
_set_prop = """ if (name == "{name}") {name} = value;\n"""
_get_prop = """ if (name == "{name}") return {name};\n"""
log(LogLevel.TRACE, "Calling dijitso just-in-time (JIT) compiler for SubDomain.")
lengthx = class_data['statements'][0]
lengthy = class_data['statements'][1]
members = ""
get_props = ""
set_props = ""
for k in class_data['properties']:
members += " double " + k + ";\n"
get_props += _get_prop.format(name=k)
set_props += _set_prop.format(name=k)
classname = signature
code_c = template_code.format(lengthx=lengthx,
lengthy=lengthy,
classname=classname,
math_header=_math_header,
members=members, constructor="",
get_props=get_props,
set_props=set_props)
code_h = ""
depends = []
return code_h, code_c, depends
def compile_pbc(lengthx, lengthy, properties):
mpi_comm = properties.pop("mpi_comm", MPI.comm_world)
cpp_data = {'statements': [lengthx,lengthy], 'properties': properties,
'name': 'subdomain', 'jit_generate': jit_generate}
subdomain = compile_class(cpp_data, mpi_comm=mpi_comm)
return subdomain
class CompiledPBC(cpp.mesh.SubDomain):
def __new__(cls, lengthx, lengthy, **kwargs):
properties = kwargs
return compile_pbc(f'{lengthx:.5E}', f'{lengthy:.5E}', properties)