Assemble "Constant" issue

1

Hello,

I have a variable k which is a “Constant” defined as

k = Constant(1.)

I know that is unnecessary to do it (because I can get the float(k) ) but I observed that if a print

print(assemble(k*dx(domain=mesh)))

the output is 1.2 instead of 1.

What is the problem?

2

You need to supply a minimal example that reproduces the issue, as the following code does not:

from dolfin import *

mesh = UnitSquareMesh(10, 10)
k = Constant(1)
print(assemble(k*dx(domain=mesh)))

fenics@5f53efd23667:/root/shared$ python3 mwe.py 
Calling FFC just-in-time (JIT) compiler, this may take some time.
Calling FFC just-in-time (JIT) compiler, this may take some time.
Calling FFC just-in-time (JIT) compiler, this may take some time.
Calling FFC just-in-time (JIT) compiler, this may take some time.
Calling FFC just-in-time (JIT) compiler, this may take some time.
1.0000000000000007
1 Like
3

My code constructs an 1D mesh with 4 regions. Here is the code:

from dolfin import *
import numpy


parameters['form_compiler']['optimize']           = True
parameters['form_compiler']['cpp_optimize']       = True
parameters['form_compiler']["cpp_optimize_flags"] = '-O2 -funroll-loops'
mpi_comm_world=MPI.comm_world
PETScOptions.set("mat_mumps_icntl_14",40.0)
set_log_level(40) # DEBUG option


class createDomain():
	def __init__(self):
		self.variables()

	def variables(self):
		self.nor      = 4                                # number of regions
		self.regions  = [0., 0.9 , 1.0 , 1.13333 , 1.2 ] # regions limits
		self.elements = [10, 10 , 50 , 50]               # number of elements


def createMesh(domain):

	if ( (domain.nor+1 != len(domain.regions) ) or ( len(domain.elements)+1 != len(domain.regions) ) ):
		print(colored("<", "yellow", attrs=['bold','blink']),
		colored("Something is wrong with mesh properties.", "blue", attrs=['bold']),
		colored(">", "yellow", attrs=['bold','blink']))
		sys.exit()

	mesh   = Mesh()
	editor = MeshEditor()

	vertices = []

	for i in range(domain.nor):
		if (i==0):
			lst = numpy.linspace( domain.regions[i], domain.regions[i+1] , domain.elements[i] + 1 )
			lst.tolist()
		else:
			lst = numpy.linspace( domain.regions[i], domain.regions[i+1] , domain.elements[i] + 1 )[1:]
			lst.tolist()

		vertices.extend(lst)


	topological_dim = 1
	geometrical_dim = 1

	num_local_vertices = sum(domain.elements) + 1
	num_local_cells    = sum(domain.elements)

	num_global_vertices = num_local_vertices
	num_global_cells    = num_local_cells

	editor.open(mesh, "interval" , topological_dim, geometrical_dim)
	editor.init_vertices_global(num_local_vertices, num_global_vertices)
	editor.init_cells_global(num_local_cells, num_global_cells)


	# Add vertices
	for i, vertex in enumerate(vertices):
		editor.add_vertex(i, numpy.array([vertex], dtype='float'))
	# Add cells
	for i in range(num_local_cells):
		editor.add_cell(i, numpy.array([i, i+1], dtype='uint'))

	editor.close()

	return(mesh)


# Load mesh properties
domain = createDomain()

# Create mesh
mesh   = createMesh(domain)

k = Constant(1)
print(assemble(k*dx(domain=mesh)))
4
\int_\Omega k\, \mathrm{d}\Omega = k \left| \Omega \right|

where \left| \Omega \right| is the volume of your domain, or in this case since the mesh is 1D simply the length of your domain, which is 1.2.

2 Likes
5

You are absolutely right. Thank you!