### Finding maximum of a specific solution component

448

views

0

I have a function with two components

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.

```
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

```
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

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

Please login to add an answer/comment or follow this question.