### Finding maximum of a specific solution component

448
views
0
12 months ago by
I have a function with two components

DG1  = FiniteElement('DG', triangle, deg)
elem = MixedElement([DG1, DG1])
V    = FunctionSpace(mesh, elem)
OC  = Function(V)
​

Say "OC" has components "O" and "C". How can I find the maximum value of "C" ? Is there a way to extract the array of "C" and find its maximum ? Ideally I do not want to make a copy.
Community: FEniCS Project

1
12 months ago by
mesh = UnitSquareMesh(2,2)
DG1e = FiniteElement('DG', triangle, 1)
Me   = MixedElement([DG1e]*2)
Q    = FunctionSpace(mesh, Me)

f    = Function(Q)
e    = Expression(('pow(x[0],2) + pow(x[1],2)', 'exp(x[0] + x[1])'), element=Me)

f.interpolate(e)

u,v  = f.split(deepcopy=True)

umax = u.vector().max()
vmax = v.vector().max()

umin = u.vector().min()
vmin = v.vector().min()

print umin, umax
print vmin, vmax
This does the job but is it possible to avoid making a copy, which happens with split(deepcopy=True) ?
written 12 months ago by Praveen C
I have not found a way to do this yet without deepcopy...
written 12 months ago by pf4d
2

without creating a deepcopy the same result can be obtained with

# Dofmaps and component dofs
dofs_0 = Q.sub(0).dofmap().dofs()
dofs_1 = Q.sub(1).dofmap().dofs()

umax = f.vector()[dofs_0].max()
vmax = f.vector()[dofs_1].max()

umin = f.vector()[dofs_0].min()
vmin = f.vector()[dofs_1].min()​
written 12 months ago by Hernán Mella
This is nice but does not seem to work in parallel. Do you know what needs to be changed ? Thanks.
written 10 months ago by Praveen C
This seems to work to get min/max of first component
# Dofmaps and component dofs
dof_range = Q.sub(0).dofmap().ownership_range()
dofs_0 = Q.sub(0).dofmap().dofs() - dof_range[0]
umax = f.vector()[dofs_0].max()
umin = f.vector()[dofs_0].min()
Is this really correct ?
written 10 months ago by Praveen C
For parallel, use MPI, for example, the max dimensions of the mesh! :

xmin = MPI.min(mpi_comm_world(), mesh.coordinates()[:,0].min())
xmax = MPI.max(mpi_comm_world(), mesh.coordinates()[:,0].max())
ymin = MPI.min(mpi_comm_world(), mesh.coordinates()[:,1].min())
ymax = MPI.max(mpi_comm_world(), mesh.coordinates()[:,1].max())
​
written 10 months ago by pf4d
DG elements coefficients do not correspond with the function values like Lagrange elements do. Are you sure you are getting the max and min values correctly?
written 12 months ago by Miguel
1
12 months ago by
If you investigate the structure a bit, it is easy to see that it can be achieved with

from dolfin import *
mesh = UnitSquareMesh(2,2)
DG1e = FiniteElement('DG', triangle, 1)
Me   = MixedElement([DG1e]*2)
Q    = FunctionSpace(mesh, Me)

f    = Function(Q)
e    = Expression(('pow(x[0],2) + pow(x[1],2)', 'exp(x[0] + x[1])'), element=Me)
#e2   = Expression(('0', '1'), element=Me)
#f.interpolate(e2)
#print f.vector().array()
f.interpolate(e)
print "max is %16.16e" % f.vector().array()[0:2:f.vector().size()].max()

​
Yes, that is a neat way to do this.
written 12 months ago by pf4d
except that you need to use the syntax

f.vector().array()[start:stop:step]​

so that

mesh = UnitSquareMesh(2,2)
DG1e = FiniteElement('DG', triangle, 1)
Me   = MixedElement([DG1e]*2)
Q    = FunctionSpace(mesh, Me)

f    = Function(Q)
e    = Expression(('pow(x[0],2) + pow(x[1],2)', 'exp(x[0] + x[1])'), element=Me)

f.interpolate(e)

u_a  = f.vector().array()[0:f.vector().size():f.value_rank() + 1]
v_a  = f.vector().array()[1:f.vector().size():f.value_rank() + 1]

umax = u_a.max()
vmax = v_a.max()

umin = u_a.min()
vmin = v_a.min()

print umin, umax
print vmin, vmax

written 12 months ago by pf4d