# How to find the coordinate where a vector function is max?

Hi! I am trying to find the mesh coordinate where my solution is maximum. In the following MWE, I am solving the Poisson equation with 0 BC on a disk of radius 3. To find the coordinate where the solution attains its max, I implemented the following lines of code:

``````import numpy as np
max_index = np.argmax(u.vector()[:])
print(mesh.coordinates()[max_index,:])
``````

This gives me a coordinate that doesn’t match with what I see in Paraview.

MWE:

``````from __future__ import print_function
from fenics import *
import matplotlib.pyplot as plt

from dolfin import *
import meshio
from dolfin import Mesh, XDMFFile, File, MeshValueCollection, cpp

# Optimization options for the form compiler
parameters["form_compiler"]["cpp_optimize"] = True
ffc_options = {"optimize": True, \
"eliminate_zeros": True, \
"precompute_basis_const": True, \
"precompute_ip_const": True}

import numpy as np
import meshio
from mshr import Circle, generate_mesh
from dolfin import Mesh, File, MeshFunction, Point, BoundaryMesh, SubDomain, plot, File
import matplotlib.pyplot as plt
from dolfin import *
from mshr import *
import numpy as np

class boundary(SubDomain):
def inside(self, x, on_boundary):
return on_boundary

class disk(SubDomain):
def inside(self, x, on_boundary):
return True

C = Circle(Point(0, 0), 3)

mesh = generate_mesh(C, 100)

boundary_markers = MeshFunction("size_t", mesh, mesh.topology().dim()-1, 0)
surface_markers = MeshFunction("size_t", mesh, mesh.topology().dim(), 0)
boundary().mark(boundary_markers, 2)
disk().mark(surface_markers, 1)

ds = Measure('ds', subdomain_data=boundary_markers)
dx = Measure('dx', subdomain_data=surface_markers)
n = FacetNormal(mesh)

W1 = FunctionSpace(mesh, "Lagrange", 1)
n = FacetNormal(mesh)
G , mu = 1, 0.1
u_D=Constant(0.0)

bc = DirichletBC(W1, u_D, boundary_markers, 2)

# Define variational problem
u = TrialFunction(W1)
v = TestFunction(W1)
f = Constant(-G/mu) # f=-G/mu
L = -f*v*dx

# Compute solution
u = Function(W1)
solve(a == L, u, bc)

import numpy as np
max_index = np.argmax(u.vector()[:])
print(mesh.coordinates()[max_index,:])
max_u = max(u.vector())
print(max_u)

file = File("dispex.pvd");
file << u;
``````

Results:

``````[-2.61086432  0.7678868 ]
22.4992211222971
``````

What I see in Paraview:

So, the result my script gives doesn’t match with what I see. I am expecting the max to be somewhat near the center. But `[-2.61086432 0.7678868 ]` is far off.

Can someone please tell me what’s wrong in here?

The mesh coordinates and the coordinates of the dofs are not the same,
To get the coordinates of the degrees of freedom, you need to use

``````x= W1.tabulate_dof_coordinates()
print(x[max_index])
``````
2 Likes

Thanks. It worked. But can you explain what `x= W1.tabulate_dof_coordinates()` does?

It tabulates (in other words compute) the coordinates of the degrees of freedom (the entries in u.vector()).

Imagine that you have a second order lagrange space on a triangle. Then your degrees of freedom no longer align with the mesh vertices, and therefore you need to compute/tabulate them

Hi, how do I find different maximum values for different subdomains ? I have created a post here How to calculate two max values of vector in two different subdomains?