Error in defining expression/sub domain

1

I am trying to define a user-expression to define a domain with multiple materials (2, in this case). The domain looks something like this.
0

And the 0.05.txt file looks something like the following (a 248 x 2 array of xy-coordinates):

6.028109416615006211e-01 5.515052548117607145e-01
9.596400633462656149e-01 4.778156527123709196e-01
.....
.....
5.210369931266580767e-01 7.947431392353718760e-01

Basically, I am creating a C++ string to define an Expression that takes value 0. inside the red circles and 1. otherwise. But the code fails to compile producing the following output (a sample snapshot of the error.log).

/tmp/tmpx2m4o8hv/dolfin_subdomain_2eb89048041746e695cf85fd8a456902.cpp: In member function ā€˜virtual bool dolfin::dolfin_subdomain_2eb89048041746e695cf85fd8a456902::inside(Eigen::Ref<const Eigen::Matrix<double, -1, 1> >, bool) constā€™:
/tmp/tmpx2m4o8hv/dolfin_subdomain_2eb89048041746e695cf85fd8a456902.cpp:63:0: error: lvalue required as decrement operand
          return pow(x[0]-0.602810941662,2)+pow(x[1]-0.551505254812,2) < 0.000132629119243 || pow(x[0]-0.959640063346,2) ...............
................
................
................ || pow(x[0]-0.521036993127,2)+pow(x[1]-0.794743139235,2) < 0.000132629119243;
 
/tmp/tmpx2m4o8hv/dolfin_subdomain_2eb89048041746e695cf85fd8a456902.cpp:63:0: error: lvalue required as decrement operand

MWE:

from dolfin import *
import numpy as np 
import os 

def generateCirclerandom(fname,rad):
    arr_circles = np.loadtxt(fname)    #fname contains the coordinates of the centres of the circles 
    circ_strng = []
    for i in range(arr_circles.shape[0]):
        circle_i = arr_circles[i]
        circ_strng.append('pow(x[0]-'+str(circle_i[0])+',2)+pow(x[1]-'+str(circle_i[1])+',2) < '+str(rad**2))
    tot_string = ' || '.join(u for u in circ_strng)
    return circ_strng, tot_string

L = 1.
fo = 0.05 
nps = 120
p0 = Point(0.,0.)
p1 = Point(L,L)
msh=RectangleMesh.create(comm,[p0,p1],[400,400],CellType.Type.quadrilateral) 
r_xy = L*np.sqrt(fo/np.pi/nps)
list_string, tot_string = generateCirclerandom('0.05.txt',r_xy)

Estring  = tot_string + ' ? ' + str(0.) +' : '+str(1.)
Evals = Expression(Estring,degree=2)

# ************ Also tried the following **********************

# omega_holes = CompiledSubDomain(tot_string,mpi_comm = MPI.comm_world)
# Evals = MeshFunction('double',msh,2)
# Evals.set_all(1.)
# omega_holes.mark(Evals,0.)


with XDMFFile(MPI.comm_world,'trial.xdmf') as wfil:
    wfil.write(Evals)
2

Here is another TL;DR version of the above error

import numpy as np 
from dolfin import *
import matplotlib.pyplot as plt 
comm=MPI.comm_world

def generateRandomCircle(fname,rad):
    arr_circles = np.loadtxt(fname)
    circ_string = []
    for circle in arr_circles:
        xc, yc = circle
        circ_string.append('pow(x[0]-'+str(xc)+',2)+pow(x[1]-'+str(yc)+',2) < '+str(rad**2))
    tot_string = ' || '.join(i for i in circ_string[:150])
    return tot_string

msh = UnitSquareMesh.create(comm,[600,600],CellType.Type.triangle)
Evals = MeshFunction('double',msh,2)
Evals.set_all(1.0)
nps = 120
fo=0.05
rx = np.sqrt(fo/np.pi/nps)
omega = CompiledSubDomain(generateRandomCircle('0.05.txt',rx),mpi_comm=comm)
omega.mark(Evals,0.0)
plt.figure(figsize=(8,8))
plot(Evals)
plt.show()

In the above script, if I change

' || '.join(i for i in circ_string[:150])

to

' || '.join(i for i in circ_string[:100])

it works perfectly fine to output the first 100 circles (roughly) as follows

trial_results
trial_results.png800Ɨ800 5.81 KB

As usual, if I change it to

tot_string = ' || '.join(i for i in circ_string)

it fails to compile. Any leads appreciated

3

Hi, consider the following (tested up to 1000 points). Bonus is that when you change the centers or radius you donā€™t have to recompile.

# The idea here is to generate code with placeholders; JIT ones and 
# then only substitute values
from dolfin import CompiledSubDomain, UnitSquareMesh, MeshFunction, File
import numpy as np


def centers(npts):
    '''Names of variables for point centers'''
    pad = len(str(npts))
    X, Y = 'X%0' + str(pad) + 'd', 'Y%0' + str(pad) + 'd'
    return [X % i for i in range(npts)], [Y % i for i in range(npts)] 


def circles_code(npts):
    '''Placeholder code'''
    base = '(pow(x[0]-%s, 2) + pow(x[1]-%s, 2) < rad_2)'
    code = ' || '.join([base % xy for xy in zip(*centers(npts))])
    
    return code

npts = 1000
code = circles_code(npts)
x_center, y_center = centers(npts)
# Generate data and build substitution dictionary
x_data, y_data = np.random.rand(2, npts)
data = dict(zip(x_center, x_data))
data.update(dict(zip(y_center, y_data)))
data['rad_2'] = 0.01**2

# Build it
domain = CompiledSubDomain(code, **data)

mesh = UnitSquareMesh(500, 500)
f = MeshFunction('size_t', mesh, 2, 0)
domain.mark(f, 1)

File('domain.pvd') << f
1 Like
4

Hi, @MiroK
Thanks for the answer. I will try this out.

Any ideas what could be going wrong with the above code that I had posted ? I could not follow from reading the error.log file.
It fails to compile with both approaches (i.e. SubDomain and/or Expression).

Also, in general, trying to define circular holes (subdomains, in general) with centers lying ā€œoutsideā€ the mesh, using SubDomain should work, right ?

5

I think I understand it; a naive error in constructing the string causes it to generate C++ code to have (x[0]--<xcoord>) and x[1]--<ycoord>. Probably this is it.