Finding maximum of a specific solution component


187
views
0
4 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

2 Answers


1
4 months ago by
pf4d  
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 4 months ago by Praveen C  
I have not found a way to do this yet without deepcopy...
written 4 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 4 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 8 weeks 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 8 weeks 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 8 weeks 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 4 months ago by Miguel  
1
4 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 4 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 4 months ago by pf4d  
Please login to add an answer/comment or follow this question.

Similar posts:
Search »