[Getfem-users] Add a new FEM with polynomial and global functions

[Yassine ZAIM]

Hello dear GetFEM users,

I am trying to add a new FEM with some global functions in addition to some polynomial functions. I found this discussion about how I can do it:

https://mail.gna.org/public/getfem-users/2014-04/msg00034.html#content

However when I look  in the file "getfem_mesh_fem_global_function.cc/h" to inspire. I don't found where the basis functions are defined? as the functions:

std::stringstream s

for (int i = 0; i<size_basis; ++i) p->base()[i] =bgeot::read_base_poly(dim,s);

for the polynomial case; And also where the name of FEM is defined? as the function: add_suffix("Name", fem_element); in the polynomial case.

I tried to add my element by programming a similar class of   

template <class FUNC> class fem : public virtual_fem {...};

in which I have determined explicitly the value of base_value, grad_base_value and hess_base_value. I inherited from this class to define my element in the file getfem_fem.cc. I define the basis and DOF by the functions:

//############# code to add the basis and DOF ###############//

 base_[i] = polynomial or global basis function;//like base_value in annex

 add_node(DOF, Point);// corresponding to each basis function

I know that for this method I could program a class with the functions eval() and derivative(). But in my case I defined the functions base_value, grad_base_value and hess_base_value explicitly without need of these methods (I think). You can see my class in the annex. By this way I get a bad result.

I hope that I was clear, and I will be thankful for your help of how I can add my element correctly.

///########## Annex ##############///

  class MyFUNC : public virtual_fem {

  protected :

    std::vector<opt_long_scalar_type> base_;

  public :

    /// Gives the array of basic functions (components).

    const std::vector<opt_long_scalar_type> &base(void) const { return base_; }

    std::vector<opt_long_scalar_type> &base(void) { return base_; }

    /** Evaluates at point x, all base functions and returns the result in

        t(nb_base,target_dim) */

    void base_value(const base_node &z, base_tensor &t) const {

    //scalar_type res = 0;

      bgeot::multi_index mi(2);

      mi[1] = target_dim(); mi[0] = short_type(nb_base(0));

      t.adjust_sizes(mi);

      base_tensor::iterator it = t.begin();

      scalar_type x = *z.begin();//z[0];

      scalar_type y = *z.end();

      *it = bgeot::to_scalar(x*y);  ++it;

      *it = bgeot::to_scalar((1-x)*y);  ++it;

      *it = bgeot::to_scalar((1-x)*(1-y));  ++it;

      *it = bgeot::to_scalar(x*(1-y));  ++it;

      *it = bgeot::to_scalar((32/1281)*sqrt(2)*(1-x)*sqrt(1-x)*(12*x-66*x*x+(143/2)*x*x*x)); ++it;

      *it = bgeot::to_scalar((32/1281)*sqrt(2)*(1-y)*sqrt(1-y)*(12*y+66*y*y+(143/2)*y*y*y));

    }

    /** Evaluates at point x, the gradient of all base functions w.r.t. the

        reference element directions 0,..,dim-1 and returns the result in

        t(nb_base,target_dim,dim) */

    void grad_base_value(const base_node &z, base_tensor &t) const {

      bgeot::multi_index mi(3);

      dim_type n = dim();

      mi[2] = n; mi[1] = target_dim(); mi[0] = short_type(nb_base(0));

      t.adjust_sizes(mi);

      base_tensor::iterator it = t.begin();

      scalar_type x = *z.begin();

      scalar_type y = *z.end();

      *it = bgeot::to_scalar(y); ++it;

      *it = bgeot::to_scalar(-y); ++it;

      *it = bgeot::to_scalar(-(1-y) ); ++it;

      *it = bgeot::to_scalar((1-y)); ++it;

      *it = bgeot::to_scalar((32/1281)*sqrt(2)*(-3/2*sqrt(1-x)*(12*x-66*x*x+(143/2)*x*x*x)+(1-x)*sqrt(1-x)*(12-132*x+(429/2)*x*x))); ++it;

      *it = bgeot::to_scalar(0); ++it;

      *it = bgeot::to_scalar(x); ++it;

      *it = bgeot::to_scalar((1-x)); ++it;

      *it = bgeot::to_scalar(-(1-x)); ++it;

      *it = bgeot::to_scalar(-x); ++it;

      *it = bgeot::to_scalar(0); ++it;

      *it = bgeot::to_scalar((32/1281)*sqrt(2)*(-3/2*sqrt(1-y)*(12*y-66*y*y+(143/2)*y*y*y)+(1-y)*sqrt(1-y)*(12-132*y+(429/2)*y*y)));

    }

    /** Evaluates at point x, the hessian of all base functions w.r.t. the

        reference element directions 0,..,dim-1 and returns the result in

        t(nb_base,target_dim,dim,dim) */

    void hess_base_value(const base_node &z, base_tensor &t) const {

      bgeot::multi_index mi(4);

      dim_type n = dim();

      mi[3] = n; mi[2] = n; mi[1] = target_dim();

      mi[0] = short_type(nb_base(0));

      t.adjust_sizes(mi);

      base_tensor::iterator it = t.begin();

      scalar_type x = *z.begin();

      scalar_type y = *z.end();

      *it = bgeot::to_scalar(0); ++it;

      *it = bgeot::to_scalar(0); ++it;

      *it = bgeot::to_scalar(0); ++it;

      *it = bgeot::to_scalar(0); ++it;

      *it = bgeot::to_scalar((32/1281)*sqrt(2)*(3/4*(1/sqrt(1-x))*(12*x-66*x*x+(143/2)*x*x*x)-3*sqrt(1-x)*(12-132*x+(429/2)*x*x)+(1-x)*sqrt(1-x)*(-132+429*x))); ++it;

      *it = bgeot::to_scalar(0); ++it;

      *it = bgeot::to_scalar(1); ++it;

      *it = bgeot::to_scalar(-1); ++it;

      *it = bgeot::to_scalar(1); ++it;

      *it = bgeot::to_scalar(-1); ++it;

      *it = bgeot::to_scalar(0); ++it;

      *it = bgeot::to_scalar(0); ++it;

      *it = bgeot::to_scalar(1); ++it;

      *it = bgeot::to_scalar(-1); ++it;

      *it = bgeot::to_scalar(1); ++it;

      *it = bgeot::to_scalar(-1); ++it;

      *it = bgeot::to_scalar(0); ++it;

      *it = bgeot::to_scalar(0); ++it;

      *it = bgeot::to_scalar(0); ++it;

      *it = bgeot::to_scalar(0); ++it;

      *it = bgeot::to_scalar(0); ++it;

      *it = bgeot::to_scalar(0); ++it;

      *it = bgeot::to_scalar(0); ++it;

      *it = bgeot::to_scalar((32/1281)*sqrt(2)*(3/4*(1/sqrt(1-y))*(12*y-66*y*y+(143/2)*y*y*y)-3*sqrt(1-y)*(12-132*y+(429/2)*y*y)+(1-y)*sqrt(1-y)*(-132+429*y)));

    }

  };

--

ZAIM Yassine 

PhD Student in Applied Mathematics