getfem-commits
[Top][All Lists]
Advanced

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

[Getfem-commits] (no subject)


From: Konstantinos Poulios
Subject: [Getfem-commits] (no subject)
Date: Wed, 20 Feb 2019 07:35:46 -0500 (EST)

branch: integration-point-variables
commit 1ff720e20cba5462db15b86bbf57eb62206075ac
Author: Konstantinos Poulios <address@hidden>
Date:   Thu Dec 20 09:46:12 2018 +0100

    Enable integration point (internal) variables
---
 src/getfem/getfem_generic_assembly.h            | 21 ++++---
 src/getfem/getfem_models.h                      |  4 ++
 src/getfem_generic_assembly_compile_and_exec.cc | 69 ++++++++++++++++++---
 src/getfem_generic_assembly_semantic.cc         | 80 ++++++++++++++++++-------
 src/getfem_generic_assembly_workspace.cc        |  6 ++
 src/getfem_models.cc                            | 11 ++++
 6 files changed, 153 insertions(+), 38 deletions(-)

diff --git a/src/getfem/getfem_generic_assembly.h 
b/src/getfem/getfem_generic_assembly.h
index 5b91109..2b80c3b 100644
--- a/src/getfem/getfem_generic_assembly.h
+++ b/src/getfem/getfem_generic_assembly.h
@@ -276,7 +276,8 @@ namespace getfem {
       gmm::sub_interval I;
       const model_real_plain_vector *V;
       const im_data *imd;
-      bgeot::multi_index qdims;  // For data having a qdim != of the fem
+      bgeot::multi_index qdims;  // For data having a qdim different than
+                                 // the qdim of the fem or im_data
                                  // (dim per dof for dof data)
                                  // and for constant variables.
 
@@ -286,17 +287,18 @@ namespace getfem {
         return q;
       }
 
-      var_description(bool is_var, bool is_fem,
-                      const mesh_fem *mmf, gmm::sub_interval I_,
-                      const model_real_plain_vector *v, const im_data *imd_,
-                      size_type Q)
-        : is_variable(is_var), is_fem_dofs(is_fem), mf(mmf), I(I_), V(v),
-          imd(imd_), qdims(1) {
+      var_description(bool is_var, bool is_fem, const mesh_fem *mf_,
+                      gmm::sub_interval I_, const model_real_plain_vector *v,
+                      const im_data *imd_, size_type Q)
+        : is_variable(is_var), is_fem_dofs(is_fem), mf(mf_), I(I_), V(v),
+          imd(imd_), qdims(1)
+      {
         GMM_ASSERT1(Q > 0, "Bad dimension");
         qdims[0] = Q;
       }
       var_description() : is_variable(false), is_fem_dofs(false),
-                          mf(0), V(0), imd(0), qdims(1) { qdims[0] = 1; }
+                          mf(0), V(0), imd(0), qdims(1)
+      { qdims[0] = 1; }
     };
 
   public:
@@ -437,6 +439,9 @@ namespace getfem {
     void add_fem_variable(const std::string &name, const mesh_fem &mf,
                           const gmm::sub_interval &I,
                           const model_real_plain_vector &VV);
+    void add_im_variable(const std::string &name, const im_data &imd,
+                         const gmm::sub_interval &I,
+                         const model_real_plain_vector &VV);
     void add_fixed_size_variable(const std::string &name,
                                  const gmm::sub_interval &I,
                                  const model_real_plain_vector &VV);
diff --git a/src/getfem/getfem_models.h b/src/getfem/getfem_models.h
index b8b697c..b9acfe0 100644
--- a/src/getfem/getfem_models.h
+++ b/src/getfem/getfem_models.h
@@ -748,6 +748,10 @@ namespace getfem {
     }
 
 
+    /** Add (internal) variable, defined at integration points.*/
+    void add_im_variable(const std::string &name, const im_data &im_data,
+                         size_type niter = 1);
+
     /** Add data, defined at integration points.*/
     void add_im_data(const std::string &name, const im_data &im_data, 
size_type niter = 1);
 
diff --git a/src/getfem_generic_assembly_compile_and_exec.cc 
b/src/getfem_generic_assembly_compile_and_exec.cc
index aa189c6..62da757 100644
--- a/src/getfem_generic_assembly_compile_and_exec.cc
+++ b/src/getfem_generic_assembly_compile_and_exec.cc
@@ -4075,6 +4075,32 @@ namespace getfem {
       interpolate(interpolate_) {}
   };
 
+  struct ga_instruction_imd_vector_assembly : public ga_instruction {
+    const base_tensor &t;
+    base_vector &V;
+    const fem_interpolation_context &ctx;
+    const gmm::sub_interval &I;
+    const im_data *imd;
+    scalar_type &coeff;
+    const size_type &ipt;
+    virtual int exec() {
+      GA_DEBUG_INFO("Instruction: vector term assembly for im_data variable");
+      size_type ifirst = I.first();
+      size_type i = t.size()
+                    * imd->filtered_index_of_point(ctx.convex_num(), ipt);
+      GMM_ASSERT1(i+t.size() <= I.size(), "Internal error 
"<<i<<"+"<<t.size()<<" <= "<<I.size());
+      for (const auto &val : t.as_vector())
+        V[ifirst+(i++)] += coeff*val;
+      return 0;
+    }
+    ga_instruction_imd_vector_assembly
+    (const base_tensor &t_, base_vector &V_,
+     const fem_interpolation_context &ctx_, const gmm::sub_interval &I_,
+     const im_data *imd_, scalar_type &coeff_, const size_type &ipt_)
+    : t(t_), V(V_), ctx(ctx_), I(I_), imd(imd_), coeff(coeff_), ipt(ipt_)
+    {}
+  };
+
   struct ga_instruction_vector_assembly : public ga_instruction {
     base_tensor &t;
     base_vector &V;
@@ -4209,6 +4235,7 @@ namespace getfem {
     const gmm::sub_interval &In1, &In2;
     const mesh_fem *mfn1, *mfn2;
     const mesh_fem **mfg1, **mfg2;
+    const im_data *imd1, *imd2;
     const scalar_type &coeff, &alpha1, &alpha2;
     const size_type &nbpt, &ipt;
     base_vector elem;
@@ -4217,7 +4244,10 @@ namespace getfem {
     virtual int exec() {
       GA_DEBUG_INFO("Instruction: matrix term assembly");
       bool empty_weight = (coeff == scalar_type(0));
-      if (ipt == 0 || interpolate) {
+      bool initialize = (ipt == 0) || interpolate || imd1 || imd2;
+      bool finalize = (ipt == nbpt-1) || interpolate || imd1 || imd2;
+
+      if (initialize) {
         if (empty_weight) elem.resize(0);
         elem.resize(t.size());
         if (!empty_weight) {
@@ -4242,7 +4272,8 @@ namespace getfem {
         for (; it != ite;) *it++ += (*itt++) * e;
         // gmm::add(gmm::scaled(t.as_vector(), coeff*alpha1*alpha2), elem);
       }
-      if (ipt == nbpt-1 || interpolate) {
+
+      if (finalize) {
         const mesh_fem *pmf1 = mfg1 ? *mfg1 : mfn1;
         const mesh_fem *pmf2 = mfg2 ? *mfg2 : mfn2;
         bool reduced = (pmf1 && pmf1->is_reduced())
@@ -4260,7 +4291,13 @@ namespace getfem {
         size_type cv2 = pmf2 ? ctx2.convex_num() : s2;
         size_type N = 1;
 
-        dofs1.assign(s1, I1.first());
+        size_type ifirst1 = I1.first(), ifirst2 = I2.first();
+        if (imd1)
+          ifirst1 +=  s1*imd1->filtered_index_of_point(ctx1.convex_num(), ipt);
+        if (imd2)
+          ifirst2 +=  s2*imd2->filtered_index_of_point(ctx2.convex_num(), ipt);
+
+        dofs1.assign(s1, ifirst1);
         if (pmf1) {
           if (!(ctx1.is_convex_num_valid())) return 0;
           N = ctx1.N();
@@ -4280,16 +4317,16 @@ namespace getfem {
           for (size_type i=0; i < s1; ++i) dofs1[i] += i;
 
         if (pmf1 == pmf2 && cv1 == cv2) {
-          if (I1.first() == I2.first()) {
+          if (ifirst1 == ifirst2) {
             add_elem_matrix(K, dofs1, dofs1, dofs1_sort, elem, ninf*1E-14, N);
           } else {
             dofs2.resize(dofs1.size());
             for (size_type i = 0; i < dofs1.size(); ++i)
-              dofs2[i] =  dofs1[i] + I2.first() - I1.first();
+              dofs2[i] =  dofs1[i] + ifirst2 - ifirst1;
             add_elem_matrix(K, dofs1, dofs2, dofs1_sort, elem, ninf*1E-14, N);
           }
         } else {
-          dofs2.assign(s2, I2.first());
+          dofs2.assign(s2, ifirst2);
           if (pmf2) {
             if (!(ctx2.is_convex_num_valid())) return 0;
             N = std::max(N, ctx2.N());
@@ -4319,14 +4356,15 @@ namespace getfem {
      const fem_interpolation_context &ctx2_,
      const gmm::sub_interval &Ir1_, const gmm::sub_interval &In1_,
      const gmm::sub_interval &Ir2_, const gmm::sub_interval &In2_,
-     const mesh_fem *mfn1_, const mesh_fem **mfg1_,
-     const mesh_fem *mfn2_, const mesh_fem **mfg2_,
+     const mesh_fem *mfn1_, const mesh_fem **mfg1_, const im_data *imd1_,
+     const mesh_fem *mfn2_, const mesh_fem **mfg2_, const im_data *imd2_,
      const scalar_type &coeff_,
      const scalar_type &alpha2_, const scalar_type &alpha1_,
      const size_type &nbpt_, const size_type &ipt_, bool interpolate_)
       : t(t_), Kr(Kr_), Kn(Kn_), ctx1(ctx1_), ctx2(ctx2_),
         Ir1(Ir1_), Ir2(Ir2_), In1(In1_), In2(In2_),
         mfn1(mfn1_), mfn2(mfn2_), mfg1(mfg1_), mfg2(mfg2_),
+        imd1(imd1_), imd2(imd2_),
         coeff(coeff_), alpha1(alpha1_), alpha2(alpha2_),
         nbpt(nbpt_), ipt(ipt_), interpolate(interpolate_),
         dofs1(0), dofs2(0) {}
@@ -6798,6 +6836,8 @@ namespace getfem {
                               "weak form has to be a scalar quantity");
                   const mesh_fem *mf
                     = workspace.associated_mf(root->name_test1);
+                  const im_data *imd
+                    = workspace.associated_im_data(root->name_test1);
                   add_interval_to_gis(workspace, root->name_test1, gis);
 
                   if (mf) {
@@ -6831,6 +6871,14 @@ namespace getfem {
                       (root->tensor(), workspace.unreduced_vector(),
                        workspace.assembled_vector(), ctx, *Ir, *In, mf, mfg,
                        gis.coeff, gis.nbpt, gis.ipt, interpolate);
+                  } else if (imd) {
+                    GMM_ASSERT1(root->interpolate_name_test1.size() == 0,
+                                "Interpolate transformation on integration "
+                                "point variable");
+                    pgai = std::make_shared<ga_instruction_imd_vector_assembly>
+                      (root->tensor(), workspace.assembled_vector(), gis.ctx,
+                       workspace.interval_of_variable(root->name_test1),
+                       imd, gis.coeff, gis.ipt);
                   } else {
                     pgai = std::make_shared<ga_instruction_vector_assembly>
                       (root->tensor(), workspace.assembled_vector(),
@@ -6848,6 +6896,9 @@ namespace getfem {
                     *mf1=workspace.associated_mf(root->name_test1),
                     *mf2=workspace.associated_mf(root->name_test2),
                     **mfg1 = 0, **mfg2 = 0;
+                  const im_data
+                    *imd1 = workspace.associated_im_data(root->name_test1),
+                    *imd2 = workspace.associated_im_data(root->name_test2);
                   const std::string &intn1 = root->interpolate_name_test1,
                                     &intn2 = root->interpolate_name_test2;
                   bool secondary1 = intn1.size() &&
@@ -6939,7 +6990,7 @@ namespace getfem {
                     pgai = std::make_shared<ga_instruction_matrix_assembly<>>
                       (root->tensor(), workspace.unreduced_matrix(),
                        workspace.assembled_matrix(), ctx1, ctx2,
-                       *Ir1, *In1, *Ir2, *In2, mf1, mfg1, mf2, mfg2,
+                       *Ir1, *In1, *Ir2, *In2, mf1, mfg1, imd1, mf2, mfg2, 
imd2,
                        gis.coeff, *alpha1, *alpha2, gis.nbpt, gis.ipt,
                        interpolate);
                   }
diff --git a/src/getfem_generic_assembly_semantic.cc 
b/src/getfem_generic_assembly_semantic.cc
index 905dc62..f3ba309 100644
--- a/src/getfem_generic_assembly_semantic.cc
+++ b/src/getfem_generic_assembly_semantic.cc
@@ -471,13 +471,15 @@ namespace getfem {
     case GA_NODE_XFEM_MINUS_HESS_TEST: case GA_NODE_XFEM_MINUS_DIVERG_TEST:
       {
         const mesh_fem *mf = workspace.associated_mf(pnode->name);
+        const im_data *imd = workspace.associated_im_data(pnode->name);
         size_type t_type = pnode->test_function_type;
         if (t_type == 1) {
           pnode->name_test1 = pnode->name;
           pnode->interpolate_name_test1 = pnode->interpolate_name;
           pnode->interpolate_name_test2 = pnode->name_test2 = "";
-          pnode->qdim1 = (mf ? workspace.qdim(pnode->name)
-                          : gmm::vect_size(workspace.value(pnode->name)));
+          pnode->qdim1 = (mf || imd)
+                         ? workspace.qdim(pnode->name)
+                         : gmm::vect_size(workspace.value(pnode->name));
           if (option == 1)
             workspace.test1.insert
               (var_trans_pair(pnode->name_test1,
@@ -489,8 +491,9 @@ namespace getfem {
           pnode->interpolate_name_test1 = pnode->name_test1 = "";
           pnode->name_test2 = pnode->name;
           pnode->interpolate_name_test2 = pnode->interpolate_name;
-          pnode->qdim2 = (mf ? workspace.qdim(pnode->name)
-                          : gmm::vect_size(workspace.value(pnode->name)));
+          pnode->qdim2 = (mf || imd)
+                         ? workspace.qdim(pnode->name)
+                         : gmm::vect_size(workspace.value(pnode->name));
           if (option == 1)
             workspace.test2.insert
               (var_trans_pair(pnode->name_test2,
@@ -499,11 +502,11 @@ namespace getfem {
             ga_throw_error(pnode->expr, pnode->pos,
                            "Invalid null size of variable");
         }
-        if (!mf) {
-          size_type n = workspace.qdim(pnode->name);
-          if (!n)
-            ga_throw_error(pnode->expr, pnode->pos,
-                           "Invalid null size of variable");
+        size_type n = workspace.qdim(pnode->name);
+        if (!n)
+          ga_throw_error(pnode->expr, pnode->pos,
+                         "Invalid null size of variable");
+        if (!mf & !imd) { // global variable
           if (n == 1) {
             pnode->init_vector_tensor(1);
             pnode->tensor()[0] = scalar_type(1);
@@ -516,6 +519,19 @@ namespace getfem {
                 pnode->tensor()(i,j) = (i==j) ? scalar_type(1)
                                               : scalar_type(0);
           }
+        } else if (imd) {
+          bgeot::multi_index mii = workspace.qdims(pnode->name);
+          if (n == 1 && mii.size() <= 1) {
+            mii.resize(1);
+            mii[0] = 1;
+          } else
+            mii.insert(mii.begin(), n);
+          pnode->t.adjust_sizes(mii);
+          GMM_ASSERT1(pnode->tensor().size() == n*n, "Internal error");
+          auto itt = pnode->tensor().begin(); // set t equal to identity
+          for (size_type i = 0; i < n; ++i)
+            for (size_type j = 0; j < n; ++j)
+              *itt++ = (i == j) ? scalar_type(1) : scalar_type(0);
         }
       }
       break;
@@ -1731,8 +1747,8 @@ namespace getfem {
               workspace.test1.insert
                 (var_trans_pair(pnode->name_test1,
                                 pnode->interpolate_name_test1));
-            pnode->qdim1 = mf ? workspace.qdim(name)
-                              : gmm::vect_size(workspace.value(name));
+            pnode->qdim1 = (mf || imd) ? workspace.qdim(name)
+                                       : gmm::vect_size(workspace.value(name));
             if (!(pnode->qdim1))
               ga_throw_error(pnode->expr, pnode->pos,
                              "Invalid null size of variable");
@@ -1743,14 +1759,14 @@ namespace getfem {
               workspace.test2.insert
                 (var_trans_pair(pnode->name_test2,
                                 pnode->interpolate_name_test2));
-            pnode->qdim2 = mf ? workspace.qdim(name)
-                              : gmm::vect_size(workspace.value(name));
+            pnode->qdim2 = (mf || imd) ? workspace.qdim(name)
+                                       : gmm::vect_size(workspace.value(name));
             if (!(pnode->qdim2))
               ga_throw_error(pnode->expr, pnode->pos,
                              "Invalid null size of variable");
           }
 
-          if (!mf && (test || !imd)) {
+          if (!mf && !imd) { // global variable
             if (prefix_id)
               ga_throw_error(pnode->expr, pnode->pos, "Gradient, Hessian or "
                         "Divergence cannot be evaluated for fixed size data.");
@@ -1781,13 +1797,37 @@ namespace getfem {
                 gmm::copy(workspace.value(name), pnode->tensor().as_vector());
               }
             }
-          } else if (!test && imd) {
+          } else if (imd) { // im_data variable
+            size_type q = workspace.qdim(name);
+            bgeot::multi_index mii = workspace.qdims(name);
+
+            if (!q) ga_throw_error(pnode->expr, pnode->pos,
+                                   "Invalid null size of variable " << name);
+            if (mii.size() > 6)
+              ga_throw_error(pnode->expr, pnode->pos,
+                            "Tensor with too many dimensions. Limited to 6");
             if (prefix_id)
               ga_throw_error(pnode->expr, pnode->pos, "Gradient, Hessian or "
                              "Divergence cannot be evaluated for im data.");
-            pnode->node_type = GA_NODE_VAL;
-            pnode->t.adjust_sizes(workspace.qdims(name));
-          } else {
+
+            pnode->node_type = test ? GA_NODE_VAL_TEST : GA_NODE_VAL;
+
+            if (test) {
+              if (q == 1 && mii.size() <= 1) {
+                mii.resize(1);
+                mii[0] = 1;
+              } else
+                mii.insert(mii.begin(), q);
+            }
+            pnode->t.adjust_sizes(mii);
+            if (test) {
+              GMM_ASSERT1(pnode->tensor().size() == q*q, "Internal error");
+              auto itt = pnode->tensor().begin(); // set t equal to identity
+              for (size_type i = 0; i < q; ++i)
+                for (size_type j = 0; j < q; ++j)
+                  *itt++ = (i == j) ? scalar_type(1) : scalar_type(0);
+            }
+          } else { // mesh_fem variable
             size_type q = workspace.qdim(name);
             size_type n = mf->linked_mesh().dim();
             bgeot::multi_index mii = workspace.qdims(name);
@@ -1806,10 +1846,8 @@ namespace getfem {
               if (test && q == 1 && mii.size() <= 1) {
                 mii.resize(1);
                 mii[0] = 2;
-              } else if (test) {
+              } else if (test)
                 mii.insert(mii.begin(), 2);
-                pnode->t.adjust_sizes(mii);
-              }
               break;
             case 1: // grad
               pnode->node_type = test ? GA_NODE_GRAD_TEST : GA_NODE_GRAD;
diff --git a/src/getfem_generic_assembly_workspace.cc 
b/src/getfem_generic_assembly_workspace.cc
index efbc6d3..4785cd9 100644
--- a/src/getfem_generic_assembly_workspace.cc
+++ b/src/getfem_generic_assembly_workspace.cc
@@ -48,6 +48,12 @@ namespace getfem {
     variables[name] = var_description(true, true, &mf, I, &VV, 0, 1);
   }
 
+  void ga_workspace::add_im_variable
+  (const std::string &name, const im_data &imd,
+   const gmm::sub_interval &I, const model_real_plain_vector &VV) {
+    variables[name] = var_description(true, false, 0, I, &VV, &imd, 1);
+  }
+
   void ga_workspace::add_fixed_size_variable
   (const std::string &name,
    const gmm::sub_interval &I, const model_real_plain_vector &VV) {
diff --git a/src/getfem_models.cc b/src/getfem_models.cc
index e7c9315..18446ac 100644
--- a/src/getfem_models.cc
+++ b/src/getfem_models.cc
@@ -770,6 +770,17 @@ namespace getfem {
     gmm::copy(t.as_vector(), set_complex_variable(name));
   }
 
+  void model::add_im_variable(const std::string &name,
+                              const im_data &im_data,
+                              size_type niter) {
+    check_name_validity(name);
+    variables[name] = var_description(true, is_complex(), false, niter);
+    variables[name].pim_data = &im_data;
+    variables[name].set_size();
+    add_dependency(im_data);
+    act_size_to_be_done = true;
+  }
+
   void model::add_im_data(const std::string &name, const im_data &im_data,
                           size_type niter) {
     check_name_validity(name);



reply via email to

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