% FUNCTION [...] = gf_model_set(model M, [operation [, args]])
%
%   Modifies a model object.
%   
%
%   * gf_model_set(model M, 'clear')
%   Clear the model.
%
%   * gf_model_set(model M, 'add fem variable', string name, mesh_fem mf[, int niter])
%   Add a variable to the model linked to a mesh_fem. `name` is the variable
%   name and `niter` is the optional number of version of the data stored,
%   for time integration schemes.
%
%   * gf_model_set(model M, 'add variable', string name, int size[, int niter])
%   Add a variable to the model of constant size. `name` is the variable
%   name and `niter` is the optional number of version of the data stored,
%   for time integration schemes.
%
%   * gf_model_set(model M, 'resize variable', string name, int size)
%   Resize a  constant size variable of the model. `name` is the variable
%   name.
%
%   * gf_model_set(model M, 'add multiplier', string name, mesh_fem mf, string primalname[, int niter])
%   Add a particular variable linked to a fem being a multiplier with
%   respect to a primal variable. The dof will be filtered with the
%   ``gmm::range_basis`` function applied on the terms of the model
%   which link the multiplier and the primal variable. This in order to
%   retain only linearly independant constraints on the primal variable.
%   Optimized for boundary multipliers. `niter` is the optional number
%   of version of the data stored, for time integration schemes.
%
%   * gf_model_set(model M, 'add fem data', string name, mesh_fem mf[, int qdim[, int niter]])
%   Add a data to the model linked to a mesh_fem. `name` is the data name,
%   `qdim` is the optional dimension of the data over the mesh_fem and
%   `niter` is the optional number of version of the data stored,
%   for time integration schemes.
%
%   * gf_model_set(model M, 'add initialized fem data', string name, mesh_fem mf, vec V)
%   Add a data to the model linked to a mesh_fem. `name` is the data name.
%   The data is initiakized with `V`. The data can be a scalar or vector
%   field.
%
%   * gf_model_set(model M, 'add data', string name, int size[, int niter])
%   Add a data to the model of constant size. `name` is the data name
%   and `niter` is the optional number of version of the data stored,
%   for time integration schemes.
%
%   * gf_model_set(model M, 'add initialized data', string name, vec V)
%   Add a fixed size data to the model linked to a mesh_fem.
%   `name` is the data name and `V` is the value of the data.
%
%   * gf_model_set(model M, 'variable', string name, vec V[, int niter])
%   Set the value of a variable or data. `name` is the data name
%   and `niter` is the optional number of version of the data stored,
%   for time integration schemes.
%
%   * gf_model_set(model M, 'to variables', vec V)
%   Set the value of the variables of the model with the vector `V`.
%   Typically, the vector `V` results of the solve of the tangent
%   linear system (usefull to solve your problem with you own solver).
%
%   * ind = gf_model_set(model M, 'add Laplacian brick', mesh_im mim, string varname[, int region])
%   Add a Laplacian term to the model relatively to the variable `varname`.
%   If this is a vector valued variable, the Laplacian term is added
%   componentwise. `region` is an optional mesh region on which the term
%   is added. If it is not specified, it is added on the whole mesh. Return
%   the brick index in the model.
%
%   * ind = gf_model_set(model M, 'add generic elliptic brick', mesh_im mim, string varname, string dataname[, int region])
%   Add a generic elliptic term to the model relatively to the variable `varname`.
%   The shape of the elliptic term depends both on the variable and the data.
%   This corresponds to a term
%   :math:`-\text{div}(a\nabla u)`
%   where :math:`a` is the data and :math:`u` the variable. The data can be a scalar,
%   a matrix or an order four tensor. The variable can be vector valued or
%   not. If the data is a scalar or a matrix and the variable is vector
%   valued then the term is added componentwise. An order four tensor data
%   is allowed for vector valued variable only. The data can be constant or
%   describbed on a fem. Of course, when the data is a tensor describe on a
%   finite element method (a tensor field) the data can be a huge vector.
%   The components of the matrix/tensor have to be stored with the fortran
%   order (columnwise) in the data vector (compatibility with blas). The
%   symmetry of the given matrix/tensor is not verified (but assumed). If
%   this is a vector valued variable, the Laplacian term is added
%   componentwise. `region` is an optional mesh region on which the term is
%   added. If it is not specified, it is added on the whole mesh. Return the
%   brick index in the model.
%
%   * ind = gf_model_set(model M, 'add source term brick', mesh_im mim, string varname, string dataname[, int region[, string directdataname]])
%   Add a source term to the model relatively to the variable `varname`.
%   The source term is represented by the data `dataname` which could be
%   constant or described on a fem. `region` is an optional mesh region
%   on which the term is added. An additional optional data `directdataname`
%   can be provided. The corresponding data vector will be directly added
%   to the right hand side without assembly. Return the brick index in the
%   model.
%
%   * ind = gf_model_set(model M, 'add normal source term brick', mesh_im mim, string varname, string dataname, int region)
%   Add a source term on the variable `varname` on a boundary `region`.
%   This region should be a boundary. The source term is represented by the
%   data `dataname` which could be constant or described on a fem. A scalar
%   product with the outward normal unit vector to the boundary is performed.
%   The main aim of this brick is to represent a Neumann condition with a
%   vector data without performing the scalar product with the normal as a
%   pre-processing. Return the brick index in the model.
%
%   * ind = gf_model_set(model M, 'add Dirichlet condition with multipliers', mesh_im mim, string varname, mult_description, int region[, string dataname])
%   Add a Dirichlet condition on the variable `varname` and the mesh
%   region `region`. This region should be a boundary. The Dirichlet
%   condition is prescribed with a multiplier variable described by
%   `mult_description`. If `mult_description` is a string this is assumed
%   to be the variable name corresponding to the multiplier (which should be
%   first declared as a multiplier variable on the mesh region in the model).
%   If it is a finite element method (mesh_fem object) then a multiplier
%   variable will be added to the model and build on this finite element
%   method (it will be restricted to the mesh region `region` and eventually
%   some conflicting dofs with some other multiplier variables will be
%   suppressed). If it is an integer, then a  multiplier variable will be
%   added to the model and build on a classical finite element of degree
%   that integer. `dataname` is the optional right hand side of  the
%   Dirichlet condition. It could be constant or described on a fem; scalar
%   or vector valued, depending on the variable on which the Dirichlet
%   condition is prescribed. Return the brick index in the model.
%
%   * ind = gf_model_set(model M, 'add Dirichlet condition with penalization', mesh_im mim, string varname, scalar coeff, int region[, string dataname, mesh_fem mf_mult])
%   Add a Dirichlet condition on the variable `varname` and the mesh
%   region `region`. This region should be a boundary. The Dirichlet
%   condition is prescribed with penalization. The penalization coefficient
%   is initially `coeff` and will be added to the data of the model.
%   `dataname` is the optional right hand side of the Dirichlet condition.
%   It could be constant or described on a fem; scalar or vector valued,
%   depending on the variable on which the Dirichlet condition is prescribed.
%   `mf_mult` is an optional parameter which allows to weaken the
%   Dirichlet condition specifying a multiplier space.
%   Return the brick index in the model.
%
%   * ind = gf_model_set(model M, 'add generalized Dirichlet condition with multipliers', mesh_im mim, string varname, mult_description, int region, string dataname, string Hname)
%   Add a Dirichlet condition on the variable `varname` and the mesh
%   region `region`.  This version is for vector field.
%   It prescribes a condition :math:`Hu = r`
%   where `H` is a matrix field. The region should be a boundary. The Dirichlet
%   condition is prescribed with a multiplier variable described by
%   `mult_description`. If `mult_description` is a string this is assumed
%   to be the variable name corresponding to the multiplier (which should be
%   first declared as a multiplier variable on the mesh region in the model).
%   If it is a finite element method (mesh_fem object) then a multiplier
%   variable will be added to the model and build on this finite element
%   method (it will be restricted to the mesh region `region` and eventually
%   some conflicting dofs with some other multiplier variables will be
%   suppressed). If it is an integer, then a  multiplier variable will be
%   added to the model and build on a classical finite element of degree
%   that integer. `dataname` is the right hand side of  the
%   Dirichlet condition. It could be constant or described on a fem; scalar
%   or vector valued, depending on the variable on which the Dirichlet
%   condition is prescribed. `Hname` is the data
%   corresponding to the matrix field `H`.
%   Return the brick index in the model.
%
%   * ind = gf_model_set(model M, 'add generalized Dirichlet condition with penalization', mesh_im mim, string varname, scalar coeff, int region, string dataname, string Hname[, mesh_fem mf_mult])
%   Add a Dirichlet condition on the variable `varname` and the mesh
%   region `region`. This version is for vector field.
%   It prescribes a condition :math:`Hu = r`
%   where `H` is a matrix field.
%   The region should be a boundary. The Dirichlet
%   condition is prescribed with penalization. The penalization coefficient
%   is intially `coeff` and will be added to the data of the model.
%   `dataname` is the right hand side of the Dirichlet condition.
%   It could be constant or described on a fem; scalar or vector valued,
%   depending on the variable on which the Dirichlet condition is prescribed.
%   `Hname` is the data
%   corresponding to the matrix field `H`. It has to be a constant matrix
%   or described on a scalar fem.
%   `mf_mult` is an optional parameter which allows to weaken the
%   Dirichlet condition specifying a multiplier space.
%   Return the brick index in the model.
%
%   * gf_model_set(model M, 'change penalization coeff', int ind_brick, scalar coeff)
%   Change the penalization coefficient of a Dirichlet condition with
%   penalization brick. If the brick is not of this kind, this
%   function has an undefined behavior.
%
%   * ind = gf_model_set(model M, 'add Helmholtz brick', mesh_im mim, string varname, string dataname[, int region])
%   Add a Helmholtz term to the model relatively to the variable `varname`.
%   `dataname` should contain the wave number. `region` is an optional mesh
%   region on which the term is added. If it is not specified, it is added
%   on the whole mesh. Return the brick index in the model.
%
%   * ind = gf_model_set(model M, 'add Fourier Robin brick', mesh_im mim, string varname, string dataname, int region)
%   Add a Fourier-Robin term to the model relatively to the variable
%   `varname`. This corresponds to a weak term of the form
%   :math:`\int (qu).v`. `dataname`
%   should contain the parameter :math:`q` of
%   the Fourier-Robin condition. `region` is the mesh region on which
%   the term is added. Return the brick index in the model.
%
%   * ind = gf_model_set(model M, 'add constraint with multipliers', string varname, string multname, spmat B, vec L)
%   Add an additional explicit constraint on the variable `varname` thank to
%   a multiplier `multname` peviously added to the model (should be a fixed
%   size variable). The constraint is :math:`BU=L`
%   with `B` being a rectangular sparse matrix. It is possible to change
%   the constraint at any time whith the methods gf_model_set(model M, 'set private matrix')
%   and gf_model_set(model M, 'set private rhs'). Return the brick index in the model.
%
%   * ind = gf_model_set(model M, 'add constraint with penalization', string varname, scalar coeff, spmat B, vec L)
%   Add an additional explicit penalized constraint on the variable `varname`.
%   The constraint is :math`BU=L` with `B` being a rectangular sparse matrix.
%   Be aware that `B` should not contain a palin row, otherwise the whole
%   tangent matrix will be plain. It is possible to change the constraint
%   at any time whith the methods gf_model_set(model M, 'set private matrix')
%   and gf_model_set(model M, 'set private rhs'). The method
%   gf_model_set(model M, 'change penalization coeff') can be used. Return the brick
%   index in the model.
%
%   * ind = gf_model_set(model M, 'add explicit matrix', string varname1, string varname2, spmat B[, int issymmetric[, int iscoercive]])
%   Add a brick representing an explicit matrix to be added to the tangent
%   linear system relatively to the variables `varname1` and `varname2`.
%   The given matrix should have has many rows as the dimension of
%   `varname1` and as many columns as the dimension of `varname2`.
%   If the two variables are different and if `issymmetric` is set to 1
%   then the transpose of the matrix is also added to the tangent system
%   (default is 0). Set `iscoercive` to 1 if the term does not affect the
%   coercivity of the tangent system (default is 0). The matrix can be
%   changed by the command gf_model_set(model M, 'set private matrix'). Return the
%   brick index in the model.
%
%   * ind = gf_model_set(model M, 'add explicit rhs', string varname, vec L)
%   Add a brick representing an explicit right hand side to be added to
%   the right hand side of the tangent linear system relatively to the
%   variable `varname`. The given rhs should have the same size than the
%   dimension of `varname`. The rhs can be changed by the command
%   gf_model_set(model M, 'set private rhs'). Return the brick index in the model.
%
%   * gf_model_set(model M, 'set private matrix', int indbrick, spmat B)
%   For some specific bricks having an internal sparse matrix
%   (explicit bricks: 'constraint brick' and 'explicit matrix brick'),
%   set this matrix.
%
%   * gf_model_set(model M, 'set private rhs', int indbrick, vec B)
%   For some specific bricks having an internal right hand side vector
%   (explicit bricks: 'constraint brick' and 'explicit rhs brick'),
%   set this rhs.
%
%   * ind = gf_model_set(model M, 'add isotropic linearized elasticity brick', mesh_im mim, string varname, string dataname_lambda, string dataname_mu[, int region])
%   Add an isotropic linearized elasticity term to the model relatively to
%   the variable `varname`. `dataname_lambda` and `dataname_mu` should
%   contain the Lame coefficients. `region` is an optional mesh region
%   on which the term is added. If it is not specified, it is added
%   on the whole mesh. Return the brick index in the model.
%
%   * ind = gf_model_set(model M, 'add linear incompressibility brick', mesh_im mim, string varname, string multname_pressure[, int region[, string dataname_coeff]])
%   Add an linear incompressibility condition on `variable`. `multname_pressure`
%   is a variable which represent the pressure. Be aware that an inf-sup
%   condition between the finite element method describing the pressure and the
%   primal variable has to be satisfied. `region` is an optional mesh region on
%   which the term is added. If it is not specified, it is added on the whole mesh.
%   `dataname_coeff` is an optional penalization coefficient for nearly
%   incompressible elasticity for instance. In this case, it is the inverse
%   of the Lame coefficient :math:`\lambda`. Return the brick index in the model.
%
%   * ind = gf_model_set(model M, 'add nonlinear elasticity brick', mesh_im mim, string varname, string constitutive_law, string dataname[, int region])
%   Add a nonlinear elasticity term to the model relatively to the
%   variable `varname`. `lawname` is the constitutive law which
%   could be 'SaintVenant Kirchhoff', 'Mooney Rivlin' or 'Ciarlet Geymonat'.
%   IMPORTANT : if the variable is defined on a 2D mesh, the plane strain
%   approximation is automatically used.
%   `dataname` is a vector of parameters for the constitutive law. Its length
%   depends on the law. It could be a short vector of constant values or a
%   vector field described on a finite element method for variable
%   coefficients. `region` is an optional mesh region on which the term
%   is added. If it is not specified, it is added on the whole mesh. Return the
%   brick index in the model.
%
%   * ind = gf_model_set(model M, 'add elastoplasticity brick', mesh_im mim ,string projname, string varname, string datalambda, string datamu, string datathreshold, string datasigma[, int region])
%   Add a nonlinear elastoplastic term to the model relatively to the
%   variable `varname`, in small deformations, for an isotropic material
%   and for a quasistatic model. `projname` is the type of projection that
%   we want to use. For the moment, only the Von Mises projection is
%   computing that we could entering 'VM' or 'Von Mises'.
%   `datasigma` is the variable representing the constraints on the material.
%   Be carefull that `varname` and `datasigma` are composed of two iterates
%   for the time scheme needed for the Newton algorithm used.
%   Moreover, the finite element method on which `varname` is described
%   is an K ordered mesh_fem, the `datasigma` one have to be at least
%   an K-1 ordered mesh_fem.
%   `datalambda` and `datamu` are the Lame coefficients of the studied
%   material.
%   `datathreshold` is the plasticity threshold of the material.
%   The three last variable could be constants or described on the
%   same finite element method.
%   `region` is an optional mesh region on which the term is added.
%   If it is not specified, it is added on the whole mesh.
%   Return the brick index in the model.
%
%   * ind = gf_model_set(model M, 'add nonlinear incompressibility brick', mesh_im mim, string varname, string multname_pressure[, int region])
%   Add an nonlinear incompressibility condition on `variable` (for large
%   strain elasticity). `multname_pressure`
%   is a variable which represent the pressure. Be aware that an inf-sup
%   condition between the finite element method describing the pressure and the
%   primal variable has to be satisfied. `region` is an optional mesh region on
%   which the term is added. If it is not specified, it is added on the
%   whole mesh. Return the brick index in the model.
%
%   * ind = gf_model_set(model M, 'add bilaplacian brick', mesh_im mim, string varname, string dataname [, int region])
%   Add a bilaplacian brick on the variable
%   `varname` and on the mesh region `region`.
%   This represent a term :math:`\Delta(D \Delta u)`.
%   where :math:`D(x)` is a coefficient determined by `dataname` which
%   could be constant or described on a f.e.m. The corresponding weak form
%   is :math:`\int D(x)\Delta u(x) \Delta v(x) dx`.
%   Return the brick index in the model.
%
%   * ind = gf_model_set(model M, 'add Kirchhoff-Love plate brick', mesh_im mim, string varname, string dataname_D, string dataname_nu [, int region])
%   Add a bilaplacian brick on the variable
%   `varname` and on the mesh region `region`.
%   This represent a term :math:`\Delta(D \Delta u)` where :math:`D(x)`
%   is a the flexion modulus determined by `dataname_D`. The term is
%   integrated by part following a Kirchhoff-Love plate model
%   with `dataname_nu` the poisson ratio.
%   Return the brick index in the model.
%
%   * ind = gf_model_set(model M, 'add normal derivative source term brick', mesh_im mim, string varname, string dataname, int region)
%   Add a normal derivative source term brick
%   :math:`F = \int b.\partial_n v` on the variable `varname` and the
%   mesh region `region`.
%   
%   Update the right hand side of the linear system.
%   `dataname` represents `b` and `varname` represents `v`.
%   Return the brick index in the model.
%
%   * ind = gf_model_set(model M, 'add Kirchhoff-Love Neumann term brick', mesh_im mim, string varname, string dataname_M, string dataname_divM, int region)
%   Add a Neumann term brick for Kirchhoff-Love model
%   on the variable `varname` and the mesh region `region`.
%   `dataname_M` represents the bending moment tensor and  `dataname_divM`
%   its divergence.
%   Return the brick index in the model.
%
%   * ind = gf_model_set(model M, 'add normal derivative Dirichlet condition with multipliers', mesh_im mim, string varname, mult_description, int region [, string dataname, int R_must_be_derivated])
%   Add a Dirichlet condition on the normal derivative of the variable
%   `varname` and on the mesh region `region` (which should be a boundary.
%   The general form is
%   :math:`\int \partial_n u(x)v(x) = \int r(x)v(x) \forall v`
%   where :math:`r(x)` is
%   the right hand side for the Dirichlet condition (0 for
%   homogeneous conditions) and :math:`v` is in a space of multipliers
%   defined by `mult_description`.
%   If `mult_description` is a string this is assumed
%   to be the variable name corresponding to the multiplier (which should be
%   first declared as a multiplier variable on the mesh region in the model).
%   If it is a finite element method (mesh_fem object) then a multiplier
%   variable will be added to the model and build on this finite element
%   method (it will be restricted to the mesh region `region` and eventually
%   some conflicting dofs with some other multiplier variables will be
%   suppressed). If it is an integer, then a  multiplier variable will be
%   added to the model and build on a classical finite element of degree
%   that integer. `dataname` is an optional parameter which represents
%   the right hand side of the Dirichlet condition.
%   If `R_must_be_derivated` is set to `true` then the normal
%   derivative of `dataname` is considered.
%   Return the brick index in the model.
%
%   * ind = gf_model_set(model M, 'add normal derivative Dirichlet condition with penalization', mesh_im mim, string varname, scalar coeff, int region [, string dataname, int R_must_be_derivated])
%   Add a Dirichlet condition on the normal derivative of the variable
%   `varname` and on the mesh region `region` (which should be a boundary.
%   The general form is
%   :math:`\int \partial_n u(x)v(x) = \int r(x)v(x) \forall v`
%   where :math:`r(x)` is
%   the right hand side for the Dirichlet condition (0 for
%   homogeneous conditions).
%   The penalization coefficient
%   is initially `coeff` and will be added to the data of the model.
%   It can be changed with the command gf_model_set(model M, 'change penalization coeff').
%   `dataname` is an optional parameter which represents
%   the right hand side of the Dirichlet condition.
%   If `R_must_be_derivated` is set to `true` then the normal
%   derivative of `dataname` is considered.
%   Return the brick index in the model.
%
%   * ind = gf_model_set(model M, 'add mass brick', mesh_im mim, string varname[, string dataname_rho[, int region]])
%   Add mass term to the model relatively to the variable `varname`.
%   If specified, the data `dataname_rho` should contain the
%   density (1 if omitted). `region` is an optional mesh region on
%   which the term is added. If it is not specified, it
%   is added on the whole mesh. Return the brick index in the model.
%
%   * ind = gf_model_set(model M, 'add basic d on dt brick', mesh_im mim, string varnameU, string dataname_dt[, string dataname_rho[, int region]])
%   Add the standard discretization of a first order time derivative on
%   `varnameU`. The parameter `dataname_rho` is the density which could
%   be omitted (the defaul value is 1). This brick should be used in
%   addition to a time dispatcher for the other terms. Return the brick
%   index in the model.
%
%   * ind = gf_model_set(model M, 'add basic d2 on dt2 brick', mesh_im mim, string varnameU,  string datanameV, string dataname_dt, string dataname_alpha,[, string dataname_rho[, int region]])
%   Add the standard discretization of a second order time derivative
%   on `varnameU`. `datanameV` is a data represented on the same finite
%   element method as U which represents the time derivative of U. The
%   parameter `dataname_rho` is the density which could be omitted (the defaul
%   value is 1). This brick should be used in addition to a time dispatcher for
%   the other terms. The time derivative :math:`v` of the
%   variable :math:`u` is preferably computed as a
%   post-traitement which depends on each scheme. The parameter `dataname_alpha`
%   depends on the time integration scheme. Return the brick index in the model.
%
%   * gf_model_set(model M, 'add theta method dispatcher', ivec bricks_indices, string theta)
%   Add a theta-method time dispatcher to a list of bricks. For instance,
%   a matrix term :math:`K` will be replaced by
%   :math:`\theta K U^{n+1} + (1-\theta) K U^{n}`.
%   
%
%   * gf_model_set(model M, 'add midpoint dispatcher', ivec bricks_indices)
%   Add a midpoint time dispatcher to a list of bricks. For instance, a
%   nonlinear term :math:`K(U)` will be replaced by
%   :math:`K((U^{n+1} +  U^{n})/2)`.
%
%   * gf_model_set(model M, 'velocity update for order two theta method', string varnameU,  string datanameV, string dataname_dt, string dataname_theta)
%   Function which udpate the velocity :math:`v^{n+1}` after
%   the computation of the displacement :math:`u^{n+1}` and
%   before the next iteration. Specific for theta-method and when the velocity is
%   included in the data of the model.
%
%   * gf_model_set(model M, 'velocity update for Newmark scheme', int id2dt2_brick, string varnameU,  string datanameV, string dataname_dt, string dataname_twobeta, string dataname_alpha)
%   Function which udpate the velocity
%   :math:`v^{n+1}` after
%   the computation of the displacement
%   :math:`u^{n+1}` and
%   before the next iteration. Specific for Newmark scheme
%   and when the velocity is
%   included in the data of the model.*
%   This version inverts the mass matrix by a
%   conjugate gradient.
%
%   * gf_model_set(model M, 'disable bricks', ivec bricks_indices)
%   Disable a brick (the brick will no longer participate to the
%   building of the tangent linear system).
%
%   * gf_model_set(model M, 'unable bricks', ivec bricks_indices)
%   Unable a disabled brick.
%
%   * gf_model_set(model M, 'first iter')
%   To be executed before the first iteration of a time integration scheme.
%
%   * gf_model_set(model M, 'next iter')
%   To be executed at the end of each iteration of a time
%   integration scheme.
%
%   * ind = gf_model_set(model M, 'add basic contact brick', string varname_u, string multname_n[, string multname_t], string dataname_r, spmat BN[, spmat BT, string dataname_friction_coeff][, string dataname_gap[, string dataname_alpha[, int symmetrized]])
%   
%   Add a contact with  or without friction brick to the model.
%   If U is the vector
%   of degrees of freedom on which the unilateral constraint is applied,
%   the matrix `BN` have to be such that this constraint is defined by
%   :math:`B_N U \le 0`. A friction condition can be considered by adding the three
%   parameters `multname_t`, `BT` and `dataname_friction_coeff`. In this case,
%   the tangential displacement is :math:`B_T U` and the matrix `BT` should have as
%   many rows as `BN` multiplied by :math:`d-1` where :math:`d` is the domain dimension.
%   In this case also, `dataname_friction_coeff` is a data which represents the
%   coefficient of friction. It can be a scalar or a vector representing a
%   value on each contact condition.  The unilateral constraint is prescribed
%   thank to a multiplier
%   `multname_n` whose dimension should be equal to the number of rows of
%   `BN`. If a friction condition is added, it is prescribed with a
%   multiplier `multname_t` whose dimension should be equal to the number
%   of rows of `BT`. The augmentation parameter `r` should be chosen in
%   a range of
%   acceptabe values (see Getfem user documentation). `dataname_gap` is an
%   optional parameter representing the initial gap. It can be a single value
%   or a vector of value. `dataname_alpha` is an optional homogenization
%   parameter for the augmentation parameter
%   (see Getfem user documentation). The parameter `symmetrized` indicates
%   that the symmetry of the tangent matrix will be kept or not (except for
%   the part representing the coupling between contact and friction which
%   cannot be symmetrized).
%
%   * gf_model_set(model M, 'contact brick set BN', int indbrick, spmat BN)
%   Can be used to set the BN matrix of a basic contact/friction brick.
%
%   * gf_model_set(model M, 'contact brick set BT', int indbrick, spmat BT)
%   Can be used to set the BT matrix of a basic contact with
%   friction brick.
%
%   * ind = gf_model_set(model M, 'add contact with rigid obstacle brick',  mesh_im mim, string varname_u, string multname_n[, string multname_t], string dataname_r[, string dataname_friction_coeff], int region, string obstacle[,  int symmetrized])
%   
%   Add a contact with or without friction condition with a rigid obstacle
%   to the model. The condition is applied on the variable `varname_u`
%   on the boundary corresponding to `region`. The rigid obstacle should
%   be described with the string `obstacle` being a signed distance to
%   the obstacle. This string should be an expression where the coordinates
%   are 'x', 'y' in 2D and 'x', 'y', 'z' in 3D. For instance, if the rigid
%   obstacle correspond to :math:`z \le 0`, the corresponding signed distance will
%   be simply "z". `multname_n` should be a fixed size variable whose size is
%   the number of degrees of freedom on boundary `region`. It represent the
%   contact equivalent nodal forces. In order to add a friction condition
%   one has to add the `multname_t` and `dataname_friction_coeff` parameters.
%   `multname_t` should be a fixed size variable whose size is
%   the number of degrees of freedom on boundary `region` multiplied by :math:`d-1`
%   where :math:`d` is the domain dimension. It represent the friction equivalent
%   nodal forces.
%   The augmentation parameter `r` should be chosen in a
%   range of acceptabe values (close to the Young modulus of the elastic
%   body, see Getfem user documentation).  `dataname_friction_coeff` is
%   the friction coefficient. It could be a scalar or a vector of values
%   representing the friction coefficient on each contact node. The
%   parameter `symmetrized` indicates that the symmetry of the tangent
%   matrix will be kept or not. Basically, this brick compute the matrix BN
%   and the vectors gap and alpha and calls the basic contact brick.
%
%   * ind = gf_model_set(model M, 'add nonmatching meshes contact brick',  mesh_im mim1[, mesh_im mim2], string varname_u1[, string varname_u2], string multname_n[, string multname_t], string dataname_r[, string dataname_fr], int rg1, int rg2[, int slave1, int slave2,  int symmetrized])
%   
%   Add a contact with or without friction condition between two faces of
%   one or two elastic bodies. The condition is applied on the variable
%   `varname_u1` or the variables `varname_u1` and `varname_u2` depending
%   if a single or two distinct displacement fields are given. Integers
%   `rg1` and `rg2` represent the regions expected to come in contact with
%   each other. In the single displacement variable case the regions defined
%   in both `rg1` and `rg2` refer to the variable `varname_u1`. In the case
%   of two displacement variables, `rg1` refers to `varname_u1` and `rg2`
%   refers to `varname_u2`. `multname_n` should be a fixed size variable
%   whose size is the number of degrees of freedom on those regions among
%   the ones defined in `rg1` and `rg2` which are characterized as "slaves".
%   It represents the contact equivalent nodal normal forces. `multname_t`
%   should be a fixed size variable whose size corresponds to the size of
%   `multname_n` multiplied by qdim - 1 . It represents the contact
%   equivalent nodal tangent (frictional) forces. The augmentation parameter
%   `r` should be chosen in a range of acceptabe values (close to the Young
%   modulus of the elastic body, see Getfem user documentation). The
%   friction coefficient stored in the parameter `fr` is either a single
%   value or a vector of the same size as `multname_n`. The optional
%   parameters `slave1` and `slave2` declare if the regions defined in `rg1`
%   and `rg2` are correspondingly considered as "slaves". By default
%   `slave1` is true and `slave2` is false, i.e. `rg1` contains the slave
%   surfaces, while 'rg2' the master surfaces. Preferrably only one of
%   `slave1` and `slave2` is set to true. The parameter `symmetrized`
%   indicates that the symmetry of the tangent matrix will be kept or not.
%   Basically, this brick computes the matrices BN and BT and the vectors
%   gap and alpha and calls the basic contact brick.
%
%
function [varargout]=gf_model_set(varargin)
  if (nargout),
    [varargout{1:nargout}]=gf_matlab('model_set', varargin{:});
  else
    gf_matlab('model_set', varargin{:});
    if (exist('ans', 'var') == 1), varargout{1}=ans; end;
  end;
% autogenerated mfile;
