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

1

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
a = dot(grad(u), grad(v))*dx
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:

max_u
max_u1354×805 67.4 KB

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?

2

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
3

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

4

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

5

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?