getfem-commits
[Top][All Lists]
Advanced

[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

[Getfem-commits] (no subject)


From: Yves Renard
Subject: [Getfem-commits] (no subject)
Date: Fri, 29 Jun 2018 10:07:00 -0400 (EDT)

branch: devel-yves-double-integ
commit fd3c9d38798c4cd7e29b140061c518b1d140e277
Author: Yves Renard <address@hidden>
Date:   Fri Jun 29 16:05:48 2018 +0200

    Adding Secondary_domain keywork in weak form language
---
 src/getfem/getfem_generic_assembly.h               |  33 ++++++
 .../getfem_generic_assembly_compile_and_exec.h     |  18 ++-
 src/getfem/getfem_generic_assembly_tree.h          |  12 ++
 src/getfem/getfem_models.h                         |  36 +++++-
 src/getfem_generic_assembly_semantic.cc            | 132 +++++++++++++++++----
 src/getfem_generic_assembly_tree.cc                |  89 +++++++++++++-
 src/getfem_generic_assembly_workspace.cc           |  41 ++++++-
 7 files changed, 321 insertions(+), 40 deletions(-)

diff --git a/src/getfem/getfem_generic_assembly.h 
b/src/getfem/getfem_generic_assembly.h
index 0e604b9..1072a47 100644
--- a/src/getfem/getfem_generic_assembly.h
+++ b/src/getfem/getfem_generic_assembly.h
@@ -133,6 +133,29 @@ namespace getfem {
   pelementary_transformation;
 
   //=========================================================================
+  //  Virtual secundary_domain object.
+  //=========================================================================
+
+  class APIDECL virtual_secondary_domain {
+    const mesh_im &mim;
+    size_type region;
+
+  public:
+
+    // interface to be reconsidered
+    virtual void give_secondary_elements(const mesh &m, size_type cv,
+                                       std::set<size_type> &elt_lst) const = 0;
+    virtual void init(const ga_workspace &workspace) const = 0;
+    virtual void finalize() const = 0;
+
+    virtual_secondary_domain(const mesh_im &mim_, size_type region_)
+      : mim(mim_), region(region_) {}
+    virtual ~virtual_secondary_domain() {}
+  };
+
+  typedef std::shared_ptr<const virtual_secondary_domain> psecondary_domain;
+
+  //=========================================================================
   // Structure dealing with macros.
   //=========================================================================
 
@@ -320,6 +343,8 @@ namespace getfem {
     VAR_SET variables;
     std::map<std::string, pinterpolate_transformation> transformations;
     std::map<std::string, pelementary_transformation> elem_transformations;
+    std::map<std::string, psecondary_domain> secondary_domains;
+    
     std::vector<tree_description> trees;
 
     std::map<std::string, std::vector<std::string> > variable_groups;
@@ -498,7 +523,15 @@ namespace getfem {
     pelementary_transformation
     elementary_transformation(const std::string &name) const;
 
+    void add_secondary_domain(const std::string &name,
+                              psecondary_domain psecdom);
 
+    bool secondary_domain_exists(const std::string &name) const;
+
+    psecondary_domain secondary_domain(const std::string &name) const;
+
+
+    
     // extract terms
     std::string extract_constant_term(const mesh &m);
     std::string extract_order1_term(const std::string &varname);
diff --git a/src/getfem/getfem_generic_assembly_compile_and_exec.h 
b/src/getfem/getfem_generic_assembly_compile_and_exec.h
index 3802ca2..a994cfd 100644
--- a/src/getfem/getfem_generic_assembly_compile_and_exec.h
+++ b/src/getfem/getfem_generic_assembly_compile_and_exec.h
@@ -99,11 +99,13 @@ namespace getfem {
     fem_precomp_pool fp_pool;
     std::map<gauss_pt_corresp, bgeot::pstored_point_tab> neighbour_corresp;
 
-    struct region_mim : std::pair<const mesh_im *, const mesh_region *> {
-      const mesh_im* mim() const { return this->first; }
-      const mesh_region* region() const { return this->second; }
+    struct region_mim
+      : std::tuple<const mesh_im *, const mesh_region *, char *> {
+      const mesh_im* mim() const { return std::get<0>(*this); }
+      const mesh_region* region() const { return std::get<1>(*this); }
       region_mim(const mesh_im *mim_, const mesh_region *region_) :
-        std::pair<const mesh_im *, const mesh_region *>(mim_, region_) {}
+      std::tuple<const mesh_im *, const mesh_region *,char *>(mim_, region_, 0)
+      {}
     };
 
     std::map<std::string, const base_vector *> extended_vars;
@@ -133,6 +135,13 @@ namespace getfem {
       std::map<const mesh_fem *, pfem_precomp> pfps;
     };
 
+    struct secondary_domain_info {
+      papprox_integration pai;
+      bool has_ctx;
+      fem_interpolation_context ctx;
+      base_small_vector Normal;
+    };
+
     struct elementary_trans_info {
       base_matrix M;
       const mesh_fem *mf;
@@ -167,6 +176,7 @@ namespace getfem {
       std::set<std::string> transformations_der;
       std::map<std::string, interpolate_info> interpolate_infos;
       std::map<std::string, elementary_trans_info> elementary_trans_infos;
+      std::map<std::string, secondary_domain_info> secondary_domain_infos;
 
       // Instructions being executed at the first Gauss point after
       // a change of integration method only.
diff --git a/src/getfem/getfem_generic_assembly_tree.h 
b/src/getfem/getfem_generic_assembly_tree.h
index 18682c0..2095d21 100644
--- a/src/getfem/getfem_generic_assembly_tree.h
+++ b/src/getfem/getfem_generic_assembly_tree.h
@@ -91,6 +91,7 @@ namespace getfem {
     GA_INTERPOLATE, // 'Interpolate' operation
     GA_INTERPOLATE_FILTER, // 'Interpolate_filter' operation
     GA_ELEMENTARY,  // 'Elementary' operation (operation at the element level)
+    GA_SECONDARY_DOMAIN,  // For the integration on a product of two domains
     GA_XFEM_PLUS,   // �valuation on the + side of a level-set for 
fem_level_set
     GA_XFEM_MINUS,  // �valuation on the - side of a level-set for 
fem_level_set
     GA_PRINT,       // 'Print' Print the tensor
@@ -166,6 +167,17 @@ namespace getfem {
     GA_NODE_ELEMENTARY_GRAD_TEST,
     GA_NODE_ELEMENTARY_HESS_TEST,
     GA_NODE_ELEMENTARY_DIVERG_TEST,
+    GA_NODE_SECONDARY_DOMAIN,
+    GA_NODE_SECONDARY_DOMAIN_VAL,
+    GA_NODE_SECONDARY_DOMAIN_GRAD,
+    GA_NODE_SECONDARY_DOMAIN_HESS,
+    GA_NODE_SECONDARY_DOMAIN_DIVERG,
+    GA_NODE_SECONDARY_DOMAIN_VAL_TEST,
+    GA_NODE_SECONDARY_DOMAIN_GRAD_TEST,
+    GA_NODE_SECONDARY_DOMAIN_HESS_TEST,
+    GA_NODE_SECONDARY_DOMAIN_DIVERG_TEST,
+    GA_NODE_SECONDARY_DOMAIN_X,
+    GA_NODE_SECONDARY_DOMAIN_NORMAL,
     GA_NODE_XFEM_PLUS,
     GA_NODE_XFEM_PLUS_VAL,
     GA_NODE_XFEM_PLUS_GRAD,
diff --git a/src/getfem/getfem_models.h b/src/getfem/getfem_models.h
index 2027a00..8704f36 100644
--- a/src/getfem/getfem_models.h
+++ b/src/getfem/getfem_models.h
@@ -328,6 +328,7 @@ namespace getfem {
     dal::bit_vector valid_bricks, active_bricks;
     std::map<std::string, pinterpolate_transformation> transformations;
     std::map<std::string, pelementary_transformation> elem_transformations;
+    std::map<std::string, psecondary_domain> secondary_domains;
 
     // Structure dealing with time integration scheme
     int time_integration; // 0 : no, 1 : time step, 2 : init
@@ -1012,11 +1013,14 @@ namespace getfem {
     */
     virtual void next_iter();
 
-    /** Add a interpolate transformation to the model to be used with the
+    /** Add an interpolate transformation to the model to be used with the
         generic assembly.
     */
     void add_interpolate_transformation(const std::string &name,
                                         pinterpolate_transformation ptrans) {
+      if (secondary_domain_exists(name))
+        GMM_ASSERT1(false, "An secondary domain with the same "
+                    "name already exists");
       if (transformations.count(name) > 0)
         GMM_ASSERT1(name.compare("neighbour_elt"), "neighbour_elt is a "
                     "reserved interpolate transformation name");
@@ -1033,7 +1037,7 @@ namespace getfem {
       return it->second;
     }
 
-    /** Tests if `name` correpsonds to an interpolate transformation.
+    /** Tests if `name` corresponds to an interpolate transformation.
     */
     bool interpolate_transformation_exists(const std::string &name) const
     { return transformations.count(name) > 0; }
@@ -1057,11 +1061,37 @@ namespace getfem {
       return it->second;
     }
 
-    /** Tests if `name` correpsonds to an elementary transformation.
+    /** Tests if `name` corresponds to an elementary transformation.
     */
     bool elementary_transformation_exists(const std::string &name) const
     { return elem_transformations.count(name) > 0; }
 
+    
+    /** Add a secondary domain to the model to be used with the
+        generic assembly.
+    */
+    void add_secondary_domain(const std::string &name,
+                              psecondary_domain ptrans) {
+      if (interpolate_transformation_exists(name))
+        GMM_ASSERT1(false, "An interpolate transformation with the same "
+                    "name already exists");secondary_domains[name] = ptrans;
+    }
+
+    /** Get a pointer to the interpolate transformation `name`.
+    */
+    psecondary_domain
+    secondary_domain(const std::string &name) const {
+      auto  it = secondary_domains.find(name);
+      GMM_ASSERT1(it != secondary_domains.end(),
+                  "Inexistent transformation " << name);
+      return it->second;
+    }
+
+    /** Tests if `name` corresponds to an interpolate transformation.
+    */
+    bool secondary_domain_exists(const std::string &name) const
+    { return secondary_domains.count(name) > 0; }
+
     /** Gives the name of the variable of index `ind_var` of the brick
         of index `ind_brick`. */
     const std::string &varname_of_brick(size_type ind_brick,
diff --git a/src/getfem_generic_assembly_semantic.cc 
b/src/getfem_generic_assembly_semantic.cc
index fda44a7..30e87ac 100644
--- a/src/getfem_generic_assembly_semantic.cc
+++ b/src/getfem_generic_assembly_semantic.cc
@@ -66,6 +66,10 @@ namespace getfem {
         pnode->node_type == GA_NODE_ELEMENTARY_GRAD ||
         pnode->node_type == GA_NODE_ELEMENTARY_HESS ||
         pnode->node_type == GA_NODE_ELEMENTARY_DIVERG ||
+        pnode->node_type == GA_NODE_SECONDARY_DOMAIN_VAL ||
+        pnode->node_type == GA_NODE_SECONDARY_DOMAIN_GRAD ||
+        pnode->node_type == GA_NODE_SECONDARY_DOMAIN_HESS ||
+        pnode->node_type == GA_NODE_SECONDARY_DOMAIN_DIVERG ||
         pnode->node_type == GA_NODE_XFEM_PLUS_VAL ||
         pnode->node_type == GA_NODE_XFEM_PLUS_GRAD ||
         pnode->node_type == GA_NODE_XFEM_PLUS_HESS ||
@@ -129,6 +133,10 @@ namespace getfem {
                          pnode->node_type == GA_NODE_ELEMENTARY_GRAD ||
                          pnode->node_type == GA_NODE_ELEMENTARY_HESS ||
                          pnode->node_type == GA_NODE_ELEMENTARY_DIVERG);
+    bool secondary_node(pnode->node_type == GA_NODE_SECONDARY_DOMAIN_VAL ||
+                        pnode->node_type == GA_NODE_SECONDARY_DOMAIN_GRAD ||
+                        pnode->node_type == GA_NODE_SECONDARY_DOMAIN_HESS ||
+                        pnode->node_type == GA_NODE_SECONDARY_DOMAIN_DIVERG);
     bool xfem_node(pnode->node_type == GA_NODE_XFEM_PLUS_VAL ||
                    pnode->node_type == GA_NODE_XFEM_PLUS_GRAD ||
                    pnode->node_type == GA_NODE_XFEM_PLUS_HESS ||
@@ -143,7 +151,8 @@ namespace getfem {
        pnode->node_type == GA_NODE_INTERPOLATE_HESS_TEST ||
        pnode->node_type == GA_NODE_INTERPOLATE_DIVERG_TEST);
 
-    if ((plain_node || interpolate_node || elementary_node || xfem_node) &&
+    if ((plain_node || interpolate_node || secondary_node ||
+         elementary_node || xfem_node) &&
         (pnode->name.compare(varname) == 0 &&
          pnode->interpolate_name.compare(interpolatename) == 0)) marked = true;
 
@@ -253,6 +262,13 @@ namespace getfem {
        GMM_ASSERT1(false,
                    "Sorry, directional derivative do not work for the moment "
                    "with elementary transformations. Future work.");
+      case GA_NODE_SECONDARY_DOMAIN_VAL_TEST:
+      case GA_NODE_SECONDARY_DOMAIN_GRAD_TEST:
+      case GA_NODE_SECONDARY_DOMAIN_HESS_TEST:
+      case GA_NODE_SECONDARY_DOMAIN_DIVERG_TEST:
+       GMM_ASSERT1(false,
+                   "Sorry, directional derivative do not work for the moment "
+                   "with secondary domains. Future work.");
       case GA_NODE_XFEM_PLUS_VAL_TEST: case GA_NODE_XFEM_PLUS_GRAD_TEST:
       case GA_NODE_XFEM_PLUS_HESS_TEST: case GA_NODE_XFEM_PLUS_DIVERG_TEST:
       case GA_NODE_XFEM_MINUS_VAL_TEST: case GA_NODE_XFEM_MINUS_GRAD_TEST:
@@ -318,6 +334,12 @@ namespace getfem {
     case GA_NODE_INTERPOLATE_HESS: case GA_NODE_INTERPOLATE_DIVERG:
     case GA_NODE_INTERPOLATE_VAL_TEST: case GA_NODE_INTERPOLATE_GRAD_TEST:
     case GA_NODE_INTERPOLATE_HESS_TEST: case GA_NODE_INTERPOLATE_DIVERG_TEST:
+    case GA_NODE_SECONDARY_DOMAIN_VAL: case GA_NODE_SECONDARY_DOMAIN_GRAD:
+    case GA_NODE_SECONDARY_DOMAIN_HESS: case GA_NODE_SECONDARY_DOMAIN_DIVERG:
+    case GA_NODE_SECONDARY_DOMAIN_VAL_TEST:
+    case GA_NODE_SECONDARY_DOMAIN_GRAD_TEST:
+    case GA_NODE_SECONDARY_DOMAIN_HESS_TEST:
+    case GA_NODE_SECONDARY_DOMAIN_DIVERG_TEST:
       c += 1.33*(1.22+ga_hash_code(pnode->name))
         + 1.66*ga_hash_code(pnode->interpolate_name);
       break;
@@ -339,6 +361,7 @@ namespace getfem {
       c += 1.33*(1.22+ga_hash_code(pnode->name));
       break;
     case GA_NODE_INTERPOLATE_X: case GA_NODE_INTERPOLATE_NORMAL:
+    case GA_NODE_SECONDARY_DOMAIN_X: case GA_NODE_SECONDARY_DOMAIN_NORMAL:
       c += M_PI*1.33*ga_hash_code(pnode->interpolate_name);
       break;
     case GA_NODE_PREDEF_FUNC: case GA_NODE_SPEC_FUNC: case GA_NODE_OPERATOR:
@@ -403,7 +426,8 @@ namespace getfem {
     case GA_NODE_ELT_K:  case GA_NODE_ELT_B: case GA_NODE_NORMAL:
     case GA_NODE_RESHAPE: case GA_NODE_IND_MOVE_LAST: case GA_NODE_SWAP_IND:
     case GA_NODE_CONTRACT: case GA_NODE_INTERPOLATE_X:
-    case GA_NODE_INTERPOLATE_NORMAL:
+    case GA_NODE_INTERPOLATE_NORMAL: case GA_NODE_SECONDARY_DOMAIN_X:
+    case GA_NODE_SECONDARY_DOMAIN_NORMAL:
       pnode->test_function_type = 0; break;
 
     case GA_NODE_ALLINDICES: pnode->test_function_type = 0; break;
@@ -421,6 +445,8 @@ namespace getfem {
     case GA_NODE_INTERPOLATE_HESS: case GA_NODE_INTERPOLATE_DIVERG:
     case GA_NODE_ELEMENTARY_VAL: case GA_NODE_ELEMENTARY_GRAD:
     case GA_NODE_ELEMENTARY_HESS: case GA_NODE_ELEMENTARY_DIVERG:
+    case GA_NODE_SECONDARY_DOMAIN_VAL:  case GA_NODE_SECONDARY_DOMAIN_GRAD:
+    case GA_NODE_SECONDARY_DOMAIN_HESS: case GA_NODE_SECONDARY_DOMAIN_DIVERG:
     case GA_NODE_XFEM_PLUS_VAL: case GA_NODE_XFEM_PLUS_GRAD:
     case GA_NODE_XFEM_PLUS_HESS: case GA_NODE_XFEM_PLUS_DIVERG:
     case GA_NODE_XFEM_MINUS_VAL: case GA_NODE_XFEM_MINUS_GRAD:
@@ -434,6 +460,10 @@ namespace getfem {
     case GA_NODE_INTERPOLATE_DERIVATIVE:
     case GA_NODE_ELEMENTARY_VAL_TEST: case GA_NODE_ELEMENTARY_GRAD_TEST:
     case GA_NODE_ELEMENTARY_HESS_TEST: case GA_NODE_ELEMENTARY_DIVERG_TEST:
+    case GA_NODE_SECONDARY_DOMAIN_VAL_TEST:
+    case GA_NODE_SECONDARY_DOMAIN_GRAD_TEST:
+    case GA_NODE_SECONDARY_DOMAIN_HESS_TEST:
+    case GA_NODE_SECONDARY_DOMAIN_DIVERG_TEST:
     case GA_NODE_XFEM_PLUS_VAL_TEST: case GA_NODE_XFEM_PLUS_GRAD_TEST:
     case GA_NODE_XFEM_PLUS_HESS_TEST: case GA_NODE_XFEM_PLUS_DIVERG_TEST:
     case GA_NODE_XFEM_MINUS_VAL_TEST: case GA_NODE_XFEM_MINUS_GRAD_TEST:
@@ -488,14 +518,20 @@ namespace getfem {
       }
       break;
 
-    case GA_NODE_INTERPOLATE:
+    case GA_NODE_INTERPOLATE: case GA_NODE_SECONDARY_DOMAIN:
       if (pnode->name.compare("X") == 0) {
-        pnode->node_type = GA_NODE_INTERPOLATE_X;
+        if (pnode->node_type == GA_NODE_INTERPOLATE)
+          pnode->node_type = GA_NODE_INTERPOLATE_X;
+        else
+          pnode->node_type = GA_NODE_SECONDARY_DOMAIN_X;
         pnode->init_vector_tensor(meshdim);
         break;
       }
       if (pnode->name.compare("Normal") == 0) {
-        pnode->node_type = GA_NODE_INTERPOLATE_NORMAL;
+        if (pnode->node_type == GA_NODE_INTERPOLATE)
+          pnode->node_type = GA_NODE_INTERPOLATE_NORMAL;
+        else
+          pnode->node_type = GA_NODE_SECONDARY_DOMAIN_NORMAL;
         pnode->init_vector_tensor(meshdim);
         break;
       }
@@ -506,13 +542,16 @@ namespace getfem {
       {
         int ndt = ((pnode->node_type == GA_NODE_INTERPOLATE) ? 1 : 0)
           + ((pnode->node_type == GA_NODE_ELEMENTARY) ? 2 : 0)
-          + ((pnode->node_type == GA_NODE_XFEM_PLUS) ? 3 : 0)
-          + ((pnode->node_type == GA_NODE_XFEM_MINUS) ? 4 : 0);
+          + ((pnode->node_type == GA_NODE_SECONDARY_DOMAIN) ? 3 : 0)
+          + ((pnode->node_type == GA_NODE_XFEM_PLUS) ? 4 : 0)
+          + ((pnode->node_type == GA_NODE_XFEM_MINUS) ? 5 : 0);
         std::string op__name =
           (pnode->node_type == GA_NODE_INTERPOLATE) ? "Interpolation" : ""
           + (pnode->node_type == GA_NODE_ELEMENTARY) ?
-             "Elementary transformation" : ""
-          + (pnode->node_type == GA_NODE_XFEM_PLUS) ? "Xfem_plus" : ""
+             "Elementary_transformation" : ""
+          + (pnode->node_type == GA_NODE_SECONDARY_DOMAIN) ?
+             "Secondary_domain" : ""
+           + (pnode->node_type == GA_NODE_XFEM_PLUS) ? "Xfem_plus" : ""
           + (pnode->node_type == GA_NODE_XFEM_MINUS) ? "Xfem_minus" : "";
 
         std::string name = pnode->name;
@@ -570,16 +609,18 @@ namespace getfem {
             switch (ndt) {
             case 1: pnode->node_type = GA_NODE_INTERPOLATE_VAL; break;
             case 2: pnode->node_type = GA_NODE_ELEMENTARY_VAL; break;
-            case 3: pnode->node_type = GA_NODE_XFEM_PLUS_VAL; break;
-            case 4: pnode->node_type = GA_NODE_XFEM_MINUS_VAL; break;
+            case 3: pnode->node_type = GA_NODE_SECONDARY_DOMAIN_VAL; break;
+            case 4: pnode->node_type = GA_NODE_XFEM_PLUS_VAL; break;
+            case 5: pnode->node_type = GA_NODE_XFEM_MINUS_VAL; break;
             default: GMM_ASSERT1(false, "internal error");
             }
           } else {
             switch (ndt) {
             case 1: pnode->node_type = GA_NODE_INTERPOLATE_VAL_TEST; break;
             case 2: pnode->node_type = GA_NODE_ELEMENTARY_VAL_TEST; break;
-            case 3: pnode->node_type = GA_NODE_XFEM_PLUS_VAL_TEST; break;
-            case 4: pnode->node_type = GA_NODE_XFEM_MINUS_VAL_TEST; break;
+            case 3: pnode->node_type = GA_NODE_SECONDARY_DOMAIN_VAL_TEST; 
break;
+            case 4: pnode->node_type = GA_NODE_XFEM_PLUS_VAL_TEST; break;
+            case 5: pnode->node_type = GA_NODE_XFEM_MINUS_VAL_TEST; break;
             default: GMM_ASSERT1(false, "internal error");
             }
             if (q == 1 && mii.size() <= 1) {
@@ -594,8 +635,9 @@ namespace getfem {
             switch (ndt) {
             case 1: pnode->node_type = GA_NODE_INTERPOLATE_GRAD; break;
             case 2: pnode->node_type = GA_NODE_ELEMENTARY_GRAD; break;
-            case 3: pnode->node_type = GA_NODE_XFEM_PLUS_GRAD; break;
-            case 4: pnode->node_type = GA_NODE_XFEM_MINUS_GRAD; break;
+            case 3: pnode->node_type = GA_NODE_SECONDARY_DOMAIN_GRAD; break;
+            case 4: pnode->node_type = GA_NODE_XFEM_PLUS_GRAD; break;
+            case 5: pnode->node_type = GA_NODE_XFEM_MINUS_GRAD; break;
             default: GMM_ASSERT1(false, "internal error");
             }
             if (n > 1) {
@@ -606,8 +648,9 @@ namespace getfem {
             switch (ndt) {
             case 1: pnode->node_type = GA_NODE_INTERPOLATE_GRAD_TEST; break;
             case 2: pnode->node_type = GA_NODE_ELEMENTARY_GRAD_TEST; break;
-            case 3: pnode->node_type = GA_NODE_XFEM_PLUS_GRAD_TEST; break;
-            case 4: pnode->node_type = GA_NODE_XFEM_MINUS_GRAD_TEST; break;
+            case 3: pnode->node_type = 
GA_NODE_SECONDARY_DOMAIN_GRAD_TEST;break;
+            case 4: pnode->node_type = GA_NODE_XFEM_PLUS_GRAD_TEST; break;
+            case 5: pnode->node_type = GA_NODE_XFEM_MINUS_GRAD_TEST; break;
             default: GMM_ASSERT1(false, "internal error");
             }
             if (q == 1 && mii.size() <= 1) {
@@ -623,8 +666,9 @@ namespace getfem {
             switch (ndt) {
             case 1: pnode->node_type = GA_NODE_INTERPOLATE_HESS; break;
             case 2: pnode->node_type = GA_NODE_ELEMENTARY_HESS; break;
-            case 3: pnode->node_type = GA_NODE_XFEM_PLUS_HESS; break;
-            case 4: pnode->node_type = GA_NODE_XFEM_MINUS_HESS; break;
+            case 3: pnode->node_type = GA_NODE_SECONDARY_DOMAIN_HESS; break;
+            case 4: pnode->node_type = GA_NODE_XFEM_PLUS_HESS; break;
+            case 5: pnode->node_type = GA_NODE_XFEM_MINUS_HESS; break;
             default: GMM_ASSERT1(false, "internal error");
             }
             if (n > 1) {
@@ -635,8 +679,9 @@ namespace getfem {
             switch (ndt) {
             case 1: pnode->node_type = GA_NODE_INTERPOLATE_HESS_TEST; break;
             case 2: pnode->node_type = GA_NODE_ELEMENTARY_HESS_TEST; break;
-            case 3: pnode->node_type = GA_NODE_XFEM_PLUS_HESS_TEST; break;
-            case 4: pnode->node_type = GA_NODE_XFEM_MINUS_HESS_TEST; break;
+            case 3: pnode->node_type = 
GA_NODE_SECONDARY_DOMAIN_HESS_TEST;break;
+            case 4: pnode->node_type = GA_NODE_XFEM_PLUS_HESS_TEST; break;
+            case 5: pnode->node_type = GA_NODE_XFEM_MINUS_HESS_TEST; break;
             default: GMM_ASSERT1(false, "internal error");
             }
             if (q == 1 && mii.size() <= 1) {
@@ -656,8 +701,9 @@ namespace getfem {
             switch (ndt) {
             case 1: pnode->node_type = GA_NODE_INTERPOLATE_DIVERG; break;
             case 2: pnode->node_type = GA_NODE_ELEMENTARY_DIVERG; break;
-            case 3: pnode->node_type = GA_NODE_XFEM_PLUS_DIVERG; break;
-            case 4: pnode->node_type = GA_NODE_XFEM_MINUS_DIVERG; break;
+            case 3: pnode->node_type = GA_NODE_SECONDARY_DOMAIN_DIVERG;break;
+            case 4: pnode->node_type = GA_NODE_XFEM_PLUS_DIVERG; break;
+            case 5: pnode->node_type = GA_NODE_XFEM_MINUS_DIVERG; break;
             default: GMM_ASSERT1(false, "internal error");
             }
             mii.resize(1);
@@ -666,8 +712,9 @@ namespace getfem {
             switch (ndt) {
             case 1: pnode->node_type = GA_NODE_INTERPOLATE_DIVERG_TEST; break;
             case 2: pnode->node_type = GA_NODE_ELEMENTARY_DIVERG_TEST; break;
-            case 3: pnode->node_type = GA_NODE_XFEM_PLUS_DIVERG_TEST; break;
-            case 4: pnode->node_type = GA_NODE_XFEM_MINUS_DIVERG_TEST; break;
+            case 3: 
pnode->node_type=GA_NODE_SECONDARY_DOMAIN_DIVERG_TEST;break;
+            case 4: pnode->node_type = GA_NODE_XFEM_PLUS_DIVERG_TEST; break;
+            case 5: pnode->node_type = GA_NODE_XFEM_MINUS_DIVERG_TEST; break;
             default: GMM_ASSERT1(false, "internal error");
             }
             mii.resize(1);
@@ -690,6 +737,12 @@ namespace getfem {
             ga_throw_error(pnode->expr, pnode->pos,
                            "Unknown elementary transformation");
           }
+        } else if (ndt == 3) {
+          if (!(workspace.secondary_domain_exists
+                (pnode->interpolate_name))) {
+            ga_throw_error(pnode->expr, pnode->pos,
+                           "Unknown secondary domain");
+          }
         }
       }
       break;
@@ -2711,6 +2764,10 @@ namespace getfem {
 
     case GA_NODE_ELEMENTARY_VAL_TEST: case GA_NODE_ELEMENTARY_GRAD_TEST:
     case GA_NODE_ELEMENTARY_HESS_TEST: case GA_NODE_ELEMENTARY_DIVERG_TEST:
+    case GA_NODE_SECONDARY_DOMAIN_VAL_TEST:
+    case GA_NODE_SECONDARY_DOMAIN_GRAD_TEST:
+    case GA_NODE_SECONDARY_DOMAIN_HESS_TEST:
+    case GA_NODE_SECONDARY_DOMAIN_DIVERG_TEST:
     case GA_NODE_XFEM_PLUS_VAL_TEST: case GA_NODE_XFEM_PLUS_GRAD_TEST:
     case GA_NODE_XFEM_PLUS_HESS_TEST: case GA_NODE_XFEM_PLUS_DIVERG_TEST:
     case GA_NODE_XFEM_MINUS_VAL_TEST: case GA_NODE_XFEM_MINUS_GRAD_TEST:
@@ -2720,11 +2777,14 @@ namespace getfem {
     case GA_NODE_SWAP_IND: case GA_NODE_IND_MOVE_LAST: case GA_NODE_CONTRACT:
     case GA_NODE_ELT_SIZE: case GA_NODE_ELT_K: case GA_NODE_ELT_B:
     case GA_NODE_CONSTANT: case GA_NODE_X: case GA_NODE_NORMAL:
+    case GA_NODE_SECONDARY_DOMAIN_X: case GA_NODE_SECONDARY_DOMAIN_NORMAL:
     case GA_NODE_OPERATOR:
       is_constant = true; break;
 
     case GA_NODE_ELEMENTARY_VAL: case GA_NODE_ELEMENTARY_GRAD:
     case GA_NODE_ELEMENTARY_HESS: case GA_NODE_ELEMENTARY_DIVERG:
+    case GA_NODE_SECONDARY_DOMAIN_VAL: case GA_NODE_SECONDARY_DOMAIN_GRAD:
+    case GA_NODE_SECONDARY_DOMAIN_HESS: case GA_NODE_SECONDARY_DOMAIN_DIVERG:
     case GA_NODE_XFEM_PLUS_VAL: case GA_NODE_XFEM_PLUS_GRAD:
     case GA_NODE_XFEM_PLUS_HESS: case GA_NODE_XFEM_PLUS_DIVERG:
     case GA_NODE_XFEM_MINUS_VAL: case GA_NODE_XFEM_MINUS_GRAD:
@@ -3015,6 +3075,12 @@ namespace getfem {
         GMM_ASSERT1(false, "Do not know how to extract a "
                     "Neumann term with an interpolate transformation");
       break;
+    case GA_NODE_SECONDARY_DOMAIN_GRAD_TEST:
+    case GA_NODE_SECONDARY_DOMAIN_DIVERG_TEST:
+      if (pnode->name.compare(varname) == 0)
+        GMM_ASSERT1(false, "Do not know how to extract a "
+                    "Neumann term with an two-domain term");
+      break;
     default:
       break;
     }
@@ -3240,6 +3306,10 @@ namespace getfem {
                          interpolatename, order);
       break;
 
+    case GA_NODE_SECONDARY_DOMAIN_VAL:
+    case GA_NODE_SECONDARY_DOMAIN_GRAD:
+    case GA_NODE_SECONDARY_DOMAIN_HESS:
+    case GA_NODE_SECONDARY_DOMAIN_DIVERG:
     case GA_NODE_ELEMENTARY_VAL:
     case GA_NODE_ELEMENTARY_GRAD:
     case GA_NODE_ELEMENTARY_HESS:
@@ -3257,6 +3327,14 @@ namespace getfem {
         mi.push_back(pnode->tensor_proper_size(i));
       pnode->t.adjust_sizes(mi);
       switch(pnode->node_type) {
+      case GA_NODE_SECONDARY_DOMAIN_VAL:
+        pnode->node_type = GA_NODE_SECONDARY_DOMAIN_VAL_TEST; break;
+      case GA_NODE_SECONDARY_DOMAIN_GRAD:
+        pnode->node_type = GA_NODE_SECONDARY_DOMAIN_GRAD_TEST; break;
+      case GA_NODE_SECONDARY_DOMAIN_HESS:
+        pnode->node_type = GA_NODE_SECONDARY_DOMAIN_HESS_TEST; break;
+      case GA_NODE_SECONDARY_DOMAIN_DIVERG:
+        pnode->node_type = GA_NODE_SECONDARY_DOMAIN_DIVERG_TEST; break;
       case GA_NODE_ELEMENTARY_VAL:
         pnode->node_type = GA_NODE_ELEMENTARY_VAL_TEST; break;
       case GA_NODE_ELEMENTARY_GRAD:
@@ -3777,6 +3855,8 @@ namespace getfem {
 
     case GA_NODE_INTERPOLATE_HESS_TEST:
     case GA_NODE_INTERPOLATE_HESS:
+    case GA_NODE_SECONDARY_DOMAIN_HESS_TEST:
+    case GA_NODE_SECONDARY_DOMAIN_HESS:
       GMM_ASSERT1(false, "Sorry, cannot derive a hessian once more");
       break;
       
@@ -4718,6 +4798,8 @@ namespace getfem {
     case GA_NODE_INTERPOLATE_DERIVATIVE:
     case GA_NODE_ELEMENTARY_VAL: case GA_NODE_ELEMENTARY_GRAD:
     case GA_NODE_ELEMENTARY_HESS: case GA_NODE_ELEMENTARY_DIVERG:
+    case GA_NODE_SECONDARY_DOMAIN_VAL: case GA_NODE_SECONDARY_DOMAIN_GRAD:
+    case GA_NODE_SECONDARY_DOMAIN_HESS: case GA_NODE_SECONDARY_DOMAIN_DIVERG:
     case GA_NODE_XFEM_PLUS_VAL: case GA_NODE_XFEM_PLUS_GRAD:
     case GA_NODE_XFEM_PLUS_HESS: case GA_NODE_XFEM_PLUS_DIVERG:
     case GA_NODE_XFEM_MINUS_VAL: case GA_NODE_XFEM_MINUS_GRAD:
diff --git a/src/getfem_generic_assembly_tree.cc 
b/src/getfem_generic_assembly_tree.cc
index c0b7f99..7aaa926 100644
--- a/src/getfem_generic_assembly_tree.cc
+++ b/src/getfem_generic_assembly_tree.cc
@@ -147,6 +147,9 @@ namespace getfem {
       if (expr.compare(token_pos, token_length,
                        "Elementary_transformation") == 0)
         return GA_ELEMENTARY;
+      if (expr.compare(token_pos, token_length, "Secondary_domain") == 0 ||
+          expr.compare(token_pos, token_length, "Secondary_Domain") == 0)
+        return GA_SECONDARY_DOMAIN;
       if (expr.compare(token_pos, token_length, "Xfem_plus") == 0)
         return GA_XFEM_PLUS;
       if (expr.compare(token_pos, token_length, "Xfem_minus") == 0)
@@ -542,7 +545,8 @@ namespace getfem {
           pnode1->nbc1 != pnode2->nbc1)
         return false;
       break;
-    case GA_NODE_INTERPOLATE_X: case GA_NODE_INTERPOLATE_NORMAL:
+    case GA_NODE_INTERPOLATE_X: case GA_NODE_SECONDARY_DOMAIN_X:
+    case GA_NODE_INTERPOLATE_NORMAL: case GA_NODE_SECONDARY_DOMAIN_NORMAL:
       if (pnode1->interpolate_name.compare(pnode2->interpolate_name))
         return false;
       break;
@@ -554,6 +558,10 @@ namespace getfem {
     case GA_NODE_INTERPOLATE_HESS_TEST: case GA_NODE_INTERPOLATE_DIVERG_TEST:
     case GA_NODE_ELEMENTARY_VAL_TEST: case GA_NODE_ELEMENTARY_GRAD_TEST:
     case GA_NODE_ELEMENTARY_HESS_TEST: case GA_NODE_ELEMENTARY_DIVERG_TEST:
+    case GA_NODE_SECONDARY_DOMAIN_VAL_TEST:
+    case GA_NODE_SECONDARY_DOMAIN_GRAD_TEST:
+    case GA_NODE_SECONDARY_DOMAIN_HESS_TEST:
+    case GA_NODE_SECONDARY_DOMAIN_DIVERG_TEST:
     case GA_NODE_XFEM_PLUS_VAL_TEST: case GA_NODE_XFEM_PLUS_GRAD_TEST:
     case GA_NODE_XFEM_PLUS_HESS_TEST: case GA_NODE_XFEM_PLUS_DIVERG_TEST:
     case GA_NODE_XFEM_MINUS_VAL_TEST: case GA_NODE_XFEM_MINUS_GRAD_TEST:
@@ -596,6 +604,8 @@ namespace getfem {
     case GA_NODE_INTERPOLATE_HESS: case GA_NODE_INTERPOLATE_DIVERG:
     case GA_NODE_ELEMENTARY_VAL: case GA_NODE_ELEMENTARY_GRAD:
     case GA_NODE_ELEMENTARY_HESS: case GA_NODE_ELEMENTARY_DIVERG:
+    case GA_NODE_SECONDARY_DOMAIN_VAL: case GA_NODE_SECONDARY_DOMAIN_GRAD:
+    case GA_NODE_SECONDARY_DOMAIN_HESS:  case GA_NODE_SECONDARY_DOMAIN_DIVERG:
     case GA_NODE_XFEM_PLUS_VAL: case GA_NODE_XFEM_PLUS_GRAD:
     case GA_NODE_XFEM_PLUS_HESS: case GA_NODE_XFEM_PLUS_DIVERG:
     case GA_NODE_XFEM_MINUS_VAL: case GA_NODE_XFEM_MINUS_GRAD:
@@ -716,7 +726,7 @@ namespace getfem {
     if (!pnode) return;
     long prec = str.precision(16);
 
-    bool is_interpolate(false), is_elementary(false);
+    bool is_interpolate(false), is_elementary(false), is_secondary(false);
     bool is_xfem_plus(false), is_xfem_minus(false);
     switch(pnode->node_type) {
     case GA_NODE_INTERPOLATE:
@@ -746,6 +756,21 @@ namespace getfem {
       is_elementary = true;
       str << "Elementary_transformation(";
       break;
+    case GA_NODE_SECONDARY_DOMAIN:
+    case GA_NODE_SECONDARY_DOMAIN_X:
+    case GA_NODE_SECONDARY_DOMAIN_NORMAL:
+    case GA_NODE_SECONDARY_DOMAIN_VAL:
+    case GA_NODE_SECONDARY_DOMAIN_GRAD:
+    case GA_NODE_SECONDARY_DOMAIN_HESS:
+    case GA_NODE_SECONDARY_DOMAIN_DIVERG:
+    case GA_NODE_SECONDARY_DOMAIN_VAL_TEST:
+    case GA_NODE_SECONDARY_DOMAIN_GRAD_TEST:
+    case GA_NODE_SECONDARY_DOMAIN_HESS_TEST:
+    case GA_NODE_SECONDARY_DOMAIN_DIVERG_TEST:
+      str << "Secondary_domain(";
+      is_secondary = true;
+      break;
+
     case GA_NODE_XFEM_PLUS:
     case GA_NODE_XFEM_PLUS_VAL:
     case GA_NODE_XFEM_PLUS_GRAD:
@@ -778,11 +803,13 @@ namespace getfem {
     case GA_NODE_GRAD:
     case GA_NODE_INTERPOLATE_GRAD:
     case GA_NODE_ELEMENTARY_GRAD:
+    case GA_NODE_SECONDARY_DOMAIN_GRAD:
     case GA_NODE_XFEM_PLUS_GRAD:
     case GA_NODE_XFEM_MINUS_GRAD:
     case GA_NODE_GRAD_TEST:
     case GA_NODE_INTERPOLATE_GRAD_TEST:
     case GA_NODE_ELEMENTARY_GRAD_TEST:
+    case GA_NODE_SECONDARY_DOMAIN_GRAD_TEST:
     case GA_NODE_XFEM_PLUS_GRAD_TEST:
     case GA_NODE_XFEM_MINUS_GRAD_TEST:
       str << "Grad_";
@@ -790,23 +817,27 @@ namespace getfem {
     case GA_NODE_HESS:
     case GA_NODE_INTERPOLATE_HESS:
     case GA_NODE_ELEMENTARY_HESS:
+    case GA_NODE_SECONDARY_DOMAIN_HESS:
     case GA_NODE_XFEM_PLUS_HESS:
     case GA_NODE_XFEM_MINUS_HESS:
     case GA_NODE_HESS_TEST:
     case GA_NODE_INTERPOLATE_HESS_TEST:
     case GA_NODE_ELEMENTARY_HESS_TEST:
+    case GA_NODE_SECONDARY_DOMAIN_HESS_TEST:
     case GA_NODE_XFEM_PLUS_HESS_TEST:
     case GA_NODE_XFEM_MINUS_HESS_TEST:
       str << "Hess_";
       break;
     case GA_NODE_DIVERG:
     case GA_NODE_INTERPOLATE_DIVERG:
+    case GA_NODE_SECONDARY_DOMAIN_DIVERG:
     case GA_NODE_ELEMENTARY_DIVERG:
     case GA_NODE_XFEM_PLUS_DIVERG:
     case GA_NODE_XFEM_MINUS_DIVERG:
     case GA_NODE_DIVERG_TEST:
     case GA_NODE_INTERPOLATE_DIVERG_TEST:
     case GA_NODE_ELEMENTARY_DIVERG_TEST:
+    case GA_NODE_SECONDARY_DOMAIN_DIVERG_TEST:
     case GA_NODE_XFEM_PLUS_DIVERG_TEST:
     case GA_NODE_XFEM_MINUS_DIVERG_TEST:
       str << "Div_";
@@ -897,10 +928,10 @@ namespace getfem {
       else if (pnode->nbc1 != size_type(-1)) str << "," << pnode->nbc1;
       str << ")";
       break;
-    case GA_NODE_INTERPOLATE_X:
+    case GA_NODE_INTERPOLATE_X: case GA_NODE_SECONDARY_DOMAIN_X:
       str << "X";
       break;
-    case GA_NODE_INTERPOLATE_NORMAL:
+    case GA_NODE_INTERPOLATE_NORMAL: case GA_NODE_SECONDARY_DOMAIN_NORMAL:
       str << "Normal";
       break;
     case GA_NODE_INTERPOLATE_DERIVATIVE:
@@ -910,26 +941,31 @@ namespace getfem {
       break;
     case GA_NODE_INTERPOLATE:
     case GA_NODE_ELEMENTARY:
+    case GA_NODE_SECONDARY_DOMAIN:
     case GA_NODE_XFEM_PLUS:
     case GA_NODE_XFEM_MINUS:
     case GA_NODE_VAL:
     case GA_NODE_INTERPOLATE_VAL:
     case GA_NODE_ELEMENTARY_VAL:
+    case GA_NODE_SECONDARY_DOMAIN_VAL:
     case GA_NODE_XFEM_PLUS_VAL:
     case GA_NODE_XFEM_MINUS_VAL:
     case GA_NODE_GRAD:
     case GA_NODE_INTERPOLATE_GRAD:
+    case GA_NODE_SECONDARY_DOMAIN_GRAD:
     case GA_NODE_ELEMENTARY_GRAD:
     case GA_NODE_XFEM_PLUS_GRAD:
     case GA_NODE_XFEM_MINUS_GRAD:
     case GA_NODE_HESS:
     case GA_NODE_INTERPOLATE_HESS:
+    case GA_NODE_SECONDARY_DOMAIN_HESS:
     case GA_NODE_ELEMENTARY_HESS:
     case GA_NODE_XFEM_PLUS_HESS:
     case GA_NODE_XFEM_MINUS_HESS:
     case GA_NODE_DIVERG:
     case GA_NODE_INTERPOLATE_DIVERG:
     case GA_NODE_ELEMENTARY_DIVERG:
+    case GA_NODE_SECONDARY_DOMAIN_DIVERG:
     case GA_NODE_XFEM_PLUS_DIVERG:
     case GA_NODE_XFEM_MINUS_DIVERG:
       str << pnode->name;
@@ -937,21 +973,25 @@ namespace getfem {
     case GA_NODE_VAL_TEST:
     case GA_NODE_INTERPOLATE_VAL_TEST:
     case GA_NODE_ELEMENTARY_VAL_TEST:
+    case GA_NODE_SECONDARY_DOMAIN_VAL_TEST:
     case GA_NODE_XFEM_PLUS_VAL_TEST:
     case GA_NODE_XFEM_MINUS_VAL_TEST:
     case GA_NODE_GRAD_TEST:
     case GA_NODE_INTERPOLATE_GRAD_TEST:
     case GA_NODE_ELEMENTARY_GRAD_TEST:
+    case GA_NODE_SECONDARY_DOMAIN_GRAD_TEST:
     case GA_NODE_XFEM_PLUS_GRAD_TEST:
     case GA_NODE_XFEM_MINUS_GRAD_TEST:
     case GA_NODE_HESS_TEST:
     case GA_NODE_INTERPOLATE_HESS_TEST:
     case GA_NODE_ELEMENTARY_HESS_TEST:
+    case GA_NODE_SECONDARY_DOMAIN_HESS_TEST:
     case GA_NODE_XFEM_PLUS_HESS_TEST:
     case GA_NODE_XFEM_MINUS_HESS_TEST:
     case GA_NODE_DIVERG_TEST:
     case GA_NODE_INTERPOLATE_DIVERG_TEST:
     case GA_NODE_ELEMENTARY_DIVERG_TEST:
+    case GA_NODE_SECONDARY_DOMAIN_DIVERG_TEST:
     case GA_NODE_XFEM_PLUS_DIVERG_TEST:
     case GA_NODE_XFEM_MINUS_DIVERG_TEST:
       str << (pnode->test_function_type == 1 ? "Test_" : "Test2_")
@@ -1124,6 +1164,8 @@ namespace getfem {
       str << "," << pnode->interpolate_name << ")";
     else if (is_elementary)
       str << "," << pnode->elementary_name << ")";
+    else if (is_secondary)
+      str << "," << pnode->interpolate_name << ")";    
     else if (is_xfem_plus || is_xfem_minus)
       str << ")";
 
@@ -1258,6 +1300,8 @@ namespace getfem {
        case GA_NAME : pnode->node_type = GA_NODE_NAME; break;
        case GA_INTERPOLATE : pnode->node_type = GA_NODE_INTERPOLATE; break;
        case GA_ELEMENTARY : pnode->node_type = GA_NODE_ELEMENTARY; break;
+       case GA_SECONDARY_DOMAIN :
+          pnode->node_type = GA_NODE_SECONDARY_DOMAIN; break;
        case GA_XFEM_PLUS : pnode->node_type = GA_NODE_XFEM_PLUS; break;
        case GA_XFEM_MINUS: pnode->node_type = GA_NODE_XFEM_MINUS; break;
        default:break;
@@ -1361,6 +1405,7 @@ namespace getfem {
     if (pnode->node_type == GA_NODE_NAME ||
        pnode->node_type == GA_NODE_INTERPOLATE ||
        pnode->node_type == GA_NODE_ELEMENTARY ||
+       pnode->node_type == GA_NODE_SECONDARY_DOMAIN ||
        pnode->node_type == GA_NODE_XFEM_PLUS ||
        pnode->node_type == GA_NODE_XFEM_MINUS) {
       std::string name = pnode->name;
@@ -1374,6 +1419,8 @@ namespace getfem {
          case GA_NODE_NAME : pnode->op_type = GA_NAME; break;
          case GA_NODE_INTERPOLATE : pnode->op_type = GA_INTERPOLATE; break;
          case GA_NODE_ELEMENTARY : pnode->op_type = GA_ELEMENTARY; break;
+         case GA_NODE_SECONDARY_DOMAIN :
+            pnode->op_type = GA_SECONDARY_DOMAIN; break;
          case GA_NODE_XFEM_PLUS : pnode->op_type = GA_XFEM_PLUS; break;
          case GA_NODE_XFEM_MINUS: pnode->op_type = GA_XFEM_MINUS; break;
          default:break;
@@ -1605,6 +1652,40 @@ namespace getfem {
           }
           break;
 
+        case GA_SECONDARY_DOMAIN:
+          {
+            tree.add_scalar(scalar_type(0), token_pos, expr);
+            tree.current_node->node_type = GA_NODE_SECONDARY_DOMAIN;
+            t_type = ga_get_token(*expr, pos, token_pos, token_length);
+            if (t_type != GA_LPAR)
+              ga_throw_error(expr, pos-1,"Missing Secondary_domain 
arguments.");
+            t_type = ga_get_token(*expr, pos, token_pos, token_length);
+            if (t_type != GA_NAME)
+              ga_throw_error(expr, pos,
+                             "First argument of Secondary_domain should be a "
+                             "variable, test function, X or Normal.");
+            tree.current_node->name = std::string(&((*expr)[token_pos]),
+                                                  token_length);
+            
+            t_type = ga_get_token(*expr, pos, token_pos, token_length);
+            if (t_type != GA_COMMA)
+              ga_throw_error(expr, pos, "Bad format for Secondary_domain "
+                             "arguments.");
+            t_type = ga_get_token(*expr, pos, token_pos, token_length);
+            if (t_type != GA_NAME)
+              ga_throw_error(expr, pos,
+                             "Second argument of Secondary_domain should be a "
+                             "transformation name.");
+            tree.current_node->interpolate_name
+              = std::string(&((*expr)[token_pos]), token_length);
+            t_type = ga_get_token(*expr, pos, token_pos, token_length);
+            if (t_type != GA_RPAR)
+              ga_throw_error(expr, pos-1, "Missing a parenthesis after "
+                             "Secondary_domain arguments.");
+            state = 2;
+          }
+          break;
+
         case GA_XFEM_PLUS:
           {
             tree.add_scalar(scalar_type(0), token_pos, expr);
diff --git a/src/getfem_generic_assembly_workspace.cc 
b/src/getfem_generic_assembly_workspace.cc
index 4bcf856..193afc0 100644
--- a/src/getfem_generic_assembly_workspace.cc
+++ b/src/getfem_generic_assembly_workspace.cc
@@ -297,6 +297,9 @@ namespace getfem {
 
   void ga_workspace::add_interpolate_transformation
   (const std::string &name, pinterpolate_transformation ptrans) {
+     if (secondary_domain_exists(name))
+      GMM_ASSERT1(false, "An secondary domain with the same "
+                  "name already exists");
     if (transformations.find(name) != transformations.end())
       GMM_ASSERT1(name.compare("neighbour_elt"), "neighbour_elt is a "
                   "reserved interpolate transformation name");
@@ -313,8 +316,7 @@ namespace getfem {
 
   pinterpolate_transformation
   ga_workspace::interpolate_transformation(const std::string &name) const {
-    std::map<std::string, pinterpolate_transformation>::const_iterator
-      it = transformations.find(name);
+    auto  it = transformations.find(name);
     if (it != transformations.end()) return it->second;
     if (md && md->interpolate_transformation_exists(name))
       return md->interpolate_transformation(name);
@@ -334,8 +336,7 @@ namespace getfem {
 
   pelementary_transformation
   ga_workspace::elementary_transformation(const std::string &name) const {
-    std::map<std::string, pelementary_transformation>::const_iterator
-      it = elem_transformations.find(name);
+    auto  it = elem_transformations.find(name);
     if (it != elem_transformations.end()) return it->second;
     if (md && md->elementary_transformation_exists(name))
       return md->elementary_transformation(name);
@@ -345,6 +346,38 @@ namespace getfem {
     GMM_ASSERT1(false, "Inexistent elementary transformation " << name);
   }
 
+  void ga_workspace::add_secondary_domain(const std::string &name,
+                                          psecondary_domain psecdom) {
+    if (interpolate_transformation_exists(name))
+      GMM_ASSERT1(false, "An interpolate transformation with the same "
+                  "name already exists");
+    secondary_domains[name] = psecdom;
+  }
+
+  bool ga_workspace::secondary_domain_exists
+  (const std::string &name) const {
+    return (md && md->secondary_domain_exists(name)) ||
+      (parent_workspace &&
+       parent_workspace->secondary_domain_exists(name)) ||
+      (secondary_domains.find(name) != secondary_domains.end());
+  }
+
+  psecondary_domain
+  ga_workspace::secondary_domain(const std::string &name) const {
+    auto it = secondary_domains.find(name);
+    if (it != secondary_domains.end()) return it->second;
+    if (md && md->secondary_domain_exists(name))
+      return md->secondary_domain(name);
+    if (parent_workspace &&
+       parent_workspace->secondary_domain_exists(name))
+      return parent_workspace->secondary_domain(name);
+    GMM_ASSERT1(false, "Inexistent secondary domain " << name);
+  }
+
+
+
+  
+
   const mesh_region &
   ga_workspace::register_region(const mesh &m, const mesh_region &region) {
     if (&m == &dummy_mesh())



reply via email to

[Prev in Thread] Current Thread [Next in Thread]