Re: [Getfem-users] element search

[Yves Renard]

> Dear professor,
> I was trying to work between two different meshes, let us call them B
> and Omega.
> On B I have defined two variables, X and Y. This variables evolve in
> time and represent the evolution of the coordinates in Omega, in the
> sense that Omega is a fictitious domain much bigger than I need. Inside
> it I have my real domain (called Bt) that is defined by the evolution of
> the variables X and Y.
> In my problem I have to assign different conditions (i.e., a different
> source term) according to the fact that in Omega I am inside Bt or not.
> I have seen that given a mesh it is possible to search a point
> (search_point(base_node &pt)). But in my case it is difficult that X and
> Y falls exactly on a mesh node... I would like to know if it is possible
> to do a research that gives me the element touched by X and Y, and then
> assign to the nodes of the element the value I want. How can I do?
> Thank you for your help!
> Mattia

Dear Mattia,

The function 

void interpolation(const mesh_fem &mf_source,
                     mesh_trans_inv &mti,
                     const VECTU &U, VECTV &V, MAT &M,
                     int version, bool extrapolation = false);

line 319 of the file getfem_interpolation.h allows you (with version = 1) to 
have the matrix representing the interpolation of mf_source on the points 
inserted in the object mti. This means that if U is the vector of values on 
the dofs of mf_source, the values on the points in mti will be
V = M * U
where the vector V stores the corresponding values.

But it seems that you want to do the inverse operation (i.e. you have V and 
you need U). Of course, the matrix M will be non-invertible in most of the 
cases. However, you can try to use a pseudo-inverse, i.e. define U as
U = transpose(M) * inverse(M * transpose(M)) * V.
Of course it does not work if  M * transpose(M) is not invertible which may 
occurs. You can also try just conjugate gradient iterations on the linear 
system
transpose(M)*M*U = transpose(M)*V
if you start with U = 0, th CG algorithm will give you the solution of minimal 
norm.

But, to answer directly to your question, there is no function which identify 
a convex from a single node because this operation is ineficient. It needs to 
iterate on each element for each search (this would be dramatically costly).
In order to be efficient, this operation needs the elements to be sorted in a 
rtree for instance. This is why it is better to gather these operations. You 
can see how the interpolation functions are build in getfem_interpolation.h 
and try to make your own function. In particular, the object mti allows you 
to have the set of points lying in a particular element and the local 
coordinate in the corresponding reference element of each of these points.

-- 
Yves.

-------------------------------------------------------------------------
  Yves Renard (address@hidden)       tel : (33) 04.72.43.87.08
  Pole de Mathematiques,                       fax : (33) 04.72.43.85.29
  Institut Camille Jordan - CNRS UMR 5208
  INSA de Lyon, Universite de Lyon
  20, rue Albert Einstein
  69621 Villeurbanne Cedex, FRANCE
  http://math.univ-lyon1.fr/~renard
-------------------------------------------------------------------------