getfem-commits
[Top][All Lists]
Advanced

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

[Getfem-commits] r5472 - in /trunk/getfem/src: getfem/getfem_assembling.


From: Yves . Renard
Subject: [Getfem-commits] r5472 - in /trunk/getfem/src: getfem/getfem_assembling.h getfem_models.cc
Date: Tue, 08 Nov 2016 08:46:34 -0000

Author: renard
Date: Tue Nov  8 09:46:34 2016
New Revision: 5472

URL: http://svn.gna.org/viewcvs/getfem?rev=5472&view=rev
Log:
still replacing old assembly by new one

Modified:
    trunk/getfem/src/getfem/getfem_assembling.h
    trunk/getfem/src/getfem_models.cc

Modified: trunk/getfem/src/getfem/getfem_assembling.h
URL: 
http://svn.gna.org/viewcvs/getfem/trunk/getfem/src/getfem/getfem_assembling.h?rev=5472&r1=5471&r2=5472&view=diff
==============================================================================
--- trunk/getfem/src/getfem/getfem_assembling.h (original)
+++ trunk/getfem/src/getfem/getfem_assembling.h Tue Nov  8 09:46:34 2016
@@ -724,11 +724,11 @@
    const VECT &A, const mesh_region &rg,const char *assembly_description,
    std::complex<T>) {
     asm_real_or_complex_1_param_mat_(gmm::real_part(M),mim,mf_u,mf_data,
-                                        gmm::real_part(A),rg,
-                                        assembly_description, T());
+                                    gmm::real_part(A),rg,
+                                    assembly_description, T());
     asm_real_or_complex_1_param_mat_(gmm::imag_part(M),mim,mf_u,mf_data,
-                                        gmm::imag_part(A),rg,
-                                        assembly_description, T());
+                                    gmm::imag_part(A),rg,
+                                    assembly_description, T());
   }
 
   /*
@@ -790,11 +790,11 @@
    const VECT &A, const mesh_region &rg,const char *assembly_description,
    std::complex<T>) {
     asm_real_or_complex_1_param_vec_(gmm::real_part(M),mim,mf_u,mf_data,
-                                        gmm::real_part(A),rg,
-                                        assembly_description, T());
+                                    gmm::real_part(A),rg,
+                                    assembly_description, T());
     asm_real_or_complex_1_param_vec_(gmm::imag_part(M),mim,mf_u,mf_data,
-                                        gmm::imag_part(A),rg,
-                                        assembly_description, T());
+                                    gmm::imag_part(A),rg,
+                                    assembly_description, T());
   }
 
   /** 
@@ -810,6 +810,19 @@
       (M, mim, mf_u, &mf_data, F, rg, "(A*Test_u):Test2_u");
   }
     
+  /** 
+     generic mass matrix assembly with an additional constant parameter
+     (on the whole mesh or on the specified boundary) 
+     @ingroup asm
+   */
+  template<typename MAT, typename VECT>
+  void asm_mass_matrix_homogeneous_param
+  (MAT &M, const mesh_im &mim, const mesh_fem &mf_u,
+   const VECT &F, const mesh_region &rg = mesh_region::all_convexes()) {
+    asm_real_or_complex_1_param_mat
+      (M, mim, mf_u, 0, F, rg, "(A*Test_u):Test2_u");
+  }
+
   /** 
       source term (for both volumic sources and boundary (Neumann) sources).
       @ingroup asm
@@ -1107,10 +1120,6 @@
     asm_stiffness_matrix_for_laplacian(M, mim, mf, mf_data, A, rg);
   }
   
-  // -------- Before this : cleaned ----------
-
-
-
   /*
     assembly of a matrix with 1 parameter (real or complex)
     (the most common here for the assembly routines below)
@@ -1155,13 +1164,6 @@
                                 assembly_description, mf_mult, T());
   }
 
-  
-
-
-  
-
-
-
 
   /**
      assembly of @f$\int_\Omega A(x)\nabla u.\nabla address@hidden, where 
@f$A(x)@f$
@@ -1191,7 +1193,7 @@
   (MAT &M, const mesh_im &mim, const mesh_fem &mf, const mesh_fem &mf_data,
    const VECT &A, const mesh_region &rg = mesh_region::all_convexes()) {
     asm_real_or_complex_1_param_mat
-      (M,mim,mf,&mf_data,A,rg,
+      (M, mim, mf, &mf_data, A, rg,
        "(Reshape(A,meshdim,meshdim)*Grad_Test_u):Grad_Test2_u");
   }
 
@@ -1202,7 +1204,7 @@
   (MAT &M, const mesh_im &mim, const mesh_fem &mf,
    const VECT &A, const mesh_region &rg = mesh_region::all_convexes()) {
     asm_real_or_complex_1_param_mat
-      (M,mim,mf,0,A,rg,
+      (M, mim, mf, 0, A, rg,
        "(Reshape(A,meshdim,meshdim)*Grad_Test_u):Grad_Test2_u");
   }
 
@@ -1213,12 +1215,9 @@
   (MAT &M, const mesh_im &mim, const mesh_fem &mf,
    const mesh_fem &mf_data, const VECT &A, 
    const mesh_region &rg = mesh_region::all_convexes()) {
-    /* GMM_ASSERT1(mf_data.get_qdim() == 1,
-       "invalid data mesh fem (Qdim=1 required)");*/
-    asm_real_or_complex_1_param
-      (M,mim,mf,mf_data,A,rg, "a=data$1(mdim(#1),mdim(#1),#2);"
-       "M$1(#1,#1)+=comp(vGrad(#1).vGrad(#1).Base(#2))"
-       "(:,l,i,:,l,j,k).a(j,i,k)");
+    asm_real_or_complex_1_param_mat
+      (M, mim, mf, &mf_data, A, rg,
+       "(Grad_Test_u*(Reshape(A,meshdim,meshdim)')):Grad_Test2_u");
   }
 
   /** The same but with a constant matrix 
@@ -1227,54 +1226,39 @@
   void asm_stiffness_matrix_for_homogeneous_scalar_elliptic_componentwise
   (MAT &M, const mesh_im &mim, const mesh_fem &mf, const VECT &A, 
    const mesh_region &rg = mesh_region::all_convexes()) {
-    /* GMM_ASSERT1(mf_data.get_qdim() == 1,
-       "invalid data mesh fem (Qdim=1 required)");*/
-    asm_real_or_complex_1_param
-      (M,mim,mf,mf,A,rg, "a=data$1(mdim(#1),mdim(#1));"
-       "M$1(#1,#1)+=comp(vGrad(#1).vGrad(#1))"
-       "(:,l,i,:,l,j).a(j,i)");
+    asm_real_or_complex_1_param_mat
+      (M, mim, mf, 0, A, rg,
+       "(Grad_Test_u*(Reshape(A,meshdim,meshdim)')):Grad_Test2_u");
   }
 
 
   /**
      Assembly of @f$\int_\Omega A(x)\nabla u.\nabla address@hidden, where 
@f$A(x)@f$
-     is a NxNxNxN (symmetric positive definite) tensor defined on mf_data.
+     is a NxNxQxQ (symmetric positive definite) tensor defined on mf_data.
   */
   template<typename MAT, typename VECT> void
   asm_stiffness_matrix_for_vector_elliptic
   (MAT &M, const mesh_im &mim, const mesh_fem &mf, const mesh_fem &mf_data, 
    const VECT &A, const mesh_region &rg = mesh_region::all_convexes()) {
-    /* GMM_ASSERT1(mf_data.get_qdim() == 1,
-       "invalid data mesh fem (Qdim=1 required)");*/
-    /* 
-       M = a_{i,j,k,l}D_{i,j}(u)D_{k,l}(v)
-    */
-    asm_real_or_complex_1_param
-      (M,mim,mf,mf_data,A,rg, 
-       "a=data$1(qdim(#1),mdim(#1),qdim(#1),mdim(#1),#2);"
-       "t=comp(vGrad(#1).vGrad(#1).Base(#2));"
-       "M(#1,#1)+= t(:,i,j,:,k,l,p).a(i,j,k,l,p)");
+   /* M = a_{i,j,k,l}D_{i,j}(u)D_{k,l}(v) */
+    asm_real_or_complex_1_param_mat
+      (M, mim, mf, &mf_data, A, rg,
+       
"(Reshape(A,qdim(u),meshdim,qdim(u),meshdim):Grad_Test_u):Grad_Test2_u");
   }
 
   /**
      Assembly of @f$\int_\Omega A(x)\nabla u.\nabla address@hidden, where 
@f$A(x)@f$
-     is a NxNxNxN (symmetric positive definite) constant tensor.
+     is a NxNxQxQ (symmetric positive definite) constant tensor.
   */
   template<typename MAT, typename VECT> void
   asm_stiffness_matrix_for_homogeneous_vector_elliptic
   (MAT &M, const mesh_im &mim, const mesh_fem &mf,
    const VECT &A, const mesh_region &rg = mesh_region::all_convexes()) {
-    /* 
-       M = a_{i,j,k,l}D_{i,j}(u)D_{k,l}(v)
-    */
-    asm_real_or_complex_1_param
-      (M,mim,mf,mf,A,rg, 
-       "a=data$1(qdim(#1),mdim(#1),qdim(#1),mdim(#1));"
-       "t=comp(vGrad(#1).vGrad(#1));"
-       "M(#1,#1)+= t(:,i,j,:,k,l).a(i,j,k,l)");
-  }
-
-
+    /* M = a_{i,j,k,l}D_{i,j}(u)D_{k,l}(v) */
+    asm_real_or_complex_1_param_mat
+      (M, mim, mf, 0, A, rg,
+       
"(Reshape(A,qdim(u),meshdim,qdim(u),meshdim):Grad_Test_u):Grad_Test2_u");
+  }
 
   /** 
       assembly of the term @f$\int_\Omega Kuv - \nabla u.\nabla 
address@hidden, 
@@ -1288,64 +1272,13 @@
   void asm_Helmholtz(MAT &M, const mesh_im &mim, const mesh_fem &mf_u,
                     const mesh_fem &mf_data, const VECT &K_squared, 
                     const mesh_region &rg = mesh_region::all_convexes()) {
-    asm_Helmholtz(M, mim, mf_u, mf_data, K_squared,rg,
-                 typename gmm::linalg_traits<VECT>::value_type());
-  }
-
-  template<typename MAT, typename VECT, typename T>
-  void asm_Helmholtz(MAT &M, const mesh_im &mim, const mesh_fem &mf_u,
-                    const mesh_fem &mf_data,
-                    const VECT &K_squared, const mesh_region &rg, T) {
-    asm_Helmholtz_real(M, mim, mf_u, mf_data, K_squared, rg);
-  }
-
-  template<typename MAT, typename VECT, typename T>
-  void asm_Helmholtz(MAT &M, const mesh_im &mim, const mesh_fem &mf_u,
-                    const mesh_fem &mf_data, const VECT &K_squared,
-                    const mesh_region &rg, std::complex<T>) {
-    asm_Helmholtz_cplx(gmm::real_part(M), gmm::imag_part(M), mim, mf_u,
-                      mf_data, gmm::real_part(K_squared),
-                      gmm::imag_part(K_squared), rg);
-  }
-
-
-  template<typename MATr, typename MATi, typename VECTr, typename VECTi>  
-  void asm_Helmholtz_cplx(const MATr &Mr, const MATi &Mi, const mesh_im &mim,
-                         const mesh_fem &mf_u, const mesh_fem &mf_data,
-                         const VECTr &K_squaredr, const VECTi &K_squaredi, 
-                         const mesh_region &rg=mesh_region::all_convexes()) {
-    generic_assembly assem("Kr=data$1(#2); Ki=data$2(#2);"
-                          "m = comp(Base(#1).Base(#1).Base(#2)); "
-                          "M$1(#1,#1)+=sym(m(:,:,i).Kr(i) - "
-                          "comp(Grad(#1).Grad(#1))(:,i,:,i));"
-                          "M$2(#1,#1)+=sym(m(:,:,i).Ki(i));");
-    assem.push_mi(mim);
-    assem.push_mf(mf_u);
-    assem.push_mf(mf_data);
-    assem.push_data(K_squaredr); //gmm::real_part(K_squared));
-    assem.push_data(K_squaredi); //gmm::imag_part(K_squared));
-    assem.push_mat(const_cast<MATr&>(Mr));
-    assem.push_mat(const_cast<MATi&>(Mi));
-    assem.assembly(rg);
-  }
-
-  template<typename MAT, typename VECT>  
-  void asm_Helmholtz_real(const MAT &M, const mesh_im &mim,
-                         const mesh_fem &mf_u, const mesh_fem &mf_data,
-                         const VECT &K_squared, 
-                         const mesh_region &rg=mesh_region::all_convexes()) {
-    generic_assembly assem("K=data$1(#2);"
-                          "m = comp(Base(#1).Base(#1).Base(#2)); "
-                          "M$1(#1,#1)+=sym(m(:,:,i).K(i) - "
-                          "comp(Grad(#1).Grad(#1))(:,i,:,i));");
-    assem.push_mi(mim);
-    assem.push_mf(mf_u);
-    assem.push_mf(mf_data);
-    assem.push_data(K_squared);
-    assem.push_mat(const_cast<MAT&>(M));
-    assem.assembly(rg);
-  }
-
+    asm_real_or_complex_1_param_mat
+      (M, mim, mf_u, &mf_data, gmm::real_part(K_squared), rg,
+       "(A*Test_u).Test2_u - Grad_Test_u:Grad_Test2_u");
+    asm_real_or_complex_1_param_mat
+      (M, mim, mf_u, &mf_data, gmm::imag_part(K_squared), rg,
+       "(A*Test_u).Test2_u");
+  }
 
   /** 
       assembly of the term @f$\int_\Omega Kuv - \nabla u.\nabla 
address@hidden, 
@@ -1359,65 +1292,12 @@
   void asm_homogeneous_Helmholtz
   (MAT &M, const mesh_im &mim, const mesh_fem &mf_u, const VECT &K_squared, 
    const mesh_region &rg = mesh_region::all_convexes()) {
-    asm_homogeneous_Helmholtz(M, mim, mf_u, K_squared, rg,
-                 typename gmm::linalg_traits<VECT>::value_type());
-  }
-
-  template<typename MAT, typename VECT, typename T>
-  void asm_homogeneous_Helmholtz(MAT &M, const mesh_im &mim,
-                                const mesh_fem &mf_u,
-                                const VECT &K_squared,
-                                const mesh_region &rg, T) {
-    asm_homogeneous_Helmholtz_real(M, mim, mf_u, K_squared, rg);
-  }
-
-  template<typename MAT, typename VECT, typename T>
-  void asm_homogeneous_Helmholtz(MAT &M, const mesh_im &mim,
-                                const mesh_fem &mf_u,
-                                const VECT &K_squared,
-                                const mesh_region &rg, std::complex<T>) {
-    asm_homogeneous_Helmholtz_cplx(gmm::real_part(M),
-                                  gmm::imag_part(M), mim, mf_u,
-                                  gmm::real_part(K_squared),
-                                  gmm::imag_part(K_squared), rg);
-  }
-
-
-  template<typename MATr, typename MATi, typename VECTr, typename VECTi>  
-  void asm_homogeneous_Helmholtz_cplx(const MATr &Mr, const MATi &Mi,
-                                     const mesh_im &mim,
-                                     const mesh_fem &mf_u,
-                                     const VECTr &K_squaredr,
-                                     const VECTi &K_squaredi, 
-                                     const mesh_region &rg) {
-    generic_assembly assem("Kr=data$1(1); Ki=data$2(1);"
-                          "m = comp(Base(#1).Base(#1)); "
-                          "M$1(#1,#1)+=sym(m(:,:).Kr(j) - "
-                          "comp(Grad(#1).Grad(#1))(:,i,:,i));"
-                          "M$2(#1,#1)+=sym(m(:,:).Ki(j));");
-    assem.push_mi(mim);
-    assem.push_mf(mf_u);
-    assem.push_data(K_squaredr);
-    assem.push_data(K_squaredi);
-    assem.push_mat(const_cast<MATr&>(Mr));
-    assem.push_mat(const_cast<MATi&>(Mi));
-    assem.assembly(rg);
-  }
-
-  template<typename MAT, typename VECT>  
-  void asm_homogeneous_Helmholtz_real(const MAT &M, const mesh_im &mim,
-                                     const mesh_fem &mf_u,
-                                     const VECT &K_squared, 
-                                     const mesh_region &rg) {
-    generic_assembly assem("K=data(1);"
-                          "m = comp(Base(#1).Base(#1)); "
-                          "M$1(#1,#1)+=sym(m(:,:).K(j) - "
-                          "comp(Grad(#1).Grad(#1))(:,i,:,i));");
-    assem.push_mi(mim);
-    assem.push_mf(mf_u);
-    assem.push_data(K_squared);
-    assem.push_mat(const_cast<MAT&>(M));
-    assem.assembly(rg);
+    asm_real_or_complex_1_param_mat
+      (M, mim, mf_u, 0, gmm::real_part(K_squared), rg,
+       "(A*Test_u).Test2_u - Grad_Test_u:Grad_Test2_u");
+    asm_real_or_complex_1_param_mat
+      (M, mim, mf_u, 0, gmm::imag_part(K_squared), rg,
+       "(A*Test_u).Test2_u");
   }
 
   enum { ASMDIR_BUILDH = 1, ASMDIR_BUILDR = 2, ASMDIR_SIMPLIFY = 4,

Modified: trunk/getfem/src/getfem_models.cc
URL: 
http://svn.gna.org/viewcvs/getfem/trunk/getfem/src/getfem_models.cc?rev=5472&r1=5471&r2=5472&view=diff
==============================================================================
--- trunk/getfem/src/getfem_models.cc   (original)
+++ trunk/getfem/src/getfem_models.cc   Tue Nov  8 09:46:34 2016
@@ -4637,33 +4637,48 @@
         }
         GMM_TRACE2("Mass term assembly for Dirichlet condition");
         if (H_version) {
-          if (mf_H)
-            asm_real_or_complex_1_param
-              (*B, mim, mf_mult, *mf_H, *H, rg, (mf_u.get_qdim() == 1) ?
-               "F=data(#2);"
-               "M(#1,#3)+=comp(Base(#1).Base(#3).Base(#2))(:,:,i).F(i)"
-               : "F=data(qdim(#1),qdim(#1),#2);"
-               
"M(#1,#3)+=comp(vBase(#1).vBase(#3).Base(#2))(:,i,:,j,k).F(i,j,k);", &mf_u);
-          else {
-            asm_real_or_complex_1_param
-              (*B, mim, mf_mult, mf_u, *H, rg, (mf_u.get_qdim() == 1) ?
-               "F=data(1);"
-               "M(#1,#2)+=comp(Base(#1).Base(#2).F(1))"
-               : "F=data(qdim(#1),qdim(#1));"
-               "M(#1,#2)+=comp(vBase(#1).vBase(#2))(:,i,:,j).F(i,j);");
-          }
+         if (mf_u.get_qdim() == 1)
+           asm_real_or_complex_1_param_mat(*B, mim, mf_mult, mf_H, *H, rg,
+                                           "(A*Test_u).Test2_u");
+         else
+           asm_real_or_complex_1_param_mat(*B, mim, mf_mult, mf_H, *H, rg,
+                           "(Reshape(A,qdim(u),qdim(u))*Test_u).Test2_u");
+          // if (mf_H)
+          //   asm_real_or_complex_1_param
+          //     (*B, mim, mf_mult, *mf_H, *H, rg, (mf_u.get_qdim() == 1) ?
+          //      "F=data(#2);"
+          //      "M(#1,#3)+=comp(Base(#1).Base(#3).Base(#2))(:,:,i).F(i)"
+          //      : "F=data(qdim(#1),qdim(#1),#2);"
+          //      
"M(#1,#3)+=comp(vBase(#1).vBase(#3).Base(#2))(:,i,:,j,k).F(i,j,k);", &mf_u);
+          // else {
+          //   asm_real_or_complex_1_param
+          //     (*B, mim, mf_mult, mf_u, *H, rg, (mf_u.get_qdim() == 1) ?
+          //      "F=data(1);"
+          //      "M(#1,#2)+=comp(Base(#1).Base(#2).F(1))"
+          //      : "F=data(qdim(#1),qdim(#1));"
+          //      "M(#1,#2)+=comp(vBase(#1).vBase(#2))(:,i,:,j).F(i,j);");
+          // }
         }
         else if (normal_component) {
-          generic_assembly assem;
-          if (mf_mult.get_qdim() == 1)
-            assem.set("M(#2,#1)+=comp(Base(#2).vBase(#1).Normal())(:,:,i,i);");
-          else
-            
assem.set("M(#2,#1)+=comp(vBase(#2).mBase(#1).Normal())(:,i,:,i,j,j);");
-          assem.push_mi(mim);
-          assem.push_mf(mf_u);
-          assem.push_mf(mf_mult);
-          assem.push_mat(*B);
-          assem.assembly(rg);
+         ga_workspace workspace;
+         gmm::sub_interval Imult(0, mf_mult.nb_dof()), Iu(0, mf_u.nb_dof());
+         base_vector mult(mf_mult.nb_dof()), u(mf_u.nb_dof());
+         workspace.add_fem_variable("mult", mf_mult, Imult, mult);
+         workspace.add_fem_variable("u", mf_u, Iu, u);
+         workspace.add_expression("Test_mult.(Test2_u.Normal)", mim, rg);
+         workspace.set_assembled_matrix(*B);
+         workspace.assembly(2);
+
+          // generic_assembly assem;
+          // if (mf_mult.get_qdim() == 1)
+          //   
assem.set("M(#2,#1)+=comp(Base(#2).vBase(#1).Normal())(:,:,i,i);");
+          // else
+          //   
assem.set("M(#2,#1)+=comp(vBase(#2).mBase(#1).Normal())(:,i,:,i,j,j);");
+          // assem.push_mi(mim);
+          // assem.push_mf(mf_u);
+          // assem.push_mf(mf_mult);
+          // assem.push_mat(*B);
+          // assem.assembly(rg);
         } else {
           asm_mass_matrix(*B, mim, mf_mult, mf_u, rg);
         }
@@ -4796,32 +4811,49 @@
         }
         GMM_TRACE2("Mass term assembly for Dirichlet condition");
         if (H_version) {
-          if (mf_H)
-            asm_real_or_complex_1_param
-              (*B, mim, mf_mult, *mf_H, *H, rg, (mf_u.get_qdim() == 1) ?
-               "F=data(#2);"
-               "M(#1,#3)+=sym(comp(Base(#1).Base(#3).Base(#2))(:,:,i).F(i))"
-               : "F=data(qdim(#1),qdim(#1),#2);"
-               
"M(#1,#3)+=sym(comp(vBase(#1).vBase(#3).Base(#2))(:,i,:,j,k).F(i,j,k));", 
&mf_u);
-          else
-             asm_real_or_complex_1_param
-              (*B, mim, mf_mult, mf_u, *H, rg, (mf_u.get_qdim() == 1) ?
-               "F=data(1);"
-               "M(#1,#2)+=sym(comp(Base(#1).Base(#2)).F(1))"
-               : "F=data(qdim(#1),qdim(#1));"
-               "M(#1,#2)+=sym(comp(vBase(#1).vBase(#2))(:,i,:,j).F(i,j));");
+         if (mf_u.get_qdim() == 1)
+           asm_real_or_complex_1_param_mat(*B, mim, mf_mult, mf_H, *H, rg,
+                                           "(A*Test_u).Test2_u");
+         else
+           asm_real_or_complex_1_param_mat(*B, mim, mf_mult, mf_H, *H, rg,
+                           "(Reshape(A,qdim(u),qdim(u))*Test_u).Test2_u");
+          // if (mf_H)
+          //   asm_real_or_complex_1_param
+          //     (*B, mim, mf_mult, *mf_H, *H, rg, (mf_u.get_qdim() == 1) ?
+          //      "F=data(#2);"
+          //      "M(#1,#3)+=sym(comp(Base(#1).Base(#3).Base(#2))(:,:,i).F(i))"
+          //      : "F=data(qdim(#1),qdim(#1),#2);"
+          //      
"M(#1,#3)+=sym(comp(vBase(#1).vBase(#3).Base(#2))(:,i,:,j,k).F(i,j,k));", 
&mf_u);
+          // else
+          //    asm_real_or_complex_1_param
+          //     (*B, mim, mf_mult, mf_u, *H, rg, (mf_u.get_qdim() == 1) ?
+          //      "F=data(1);"
+          //      "M(#1,#2)+=sym(comp(Base(#1).Base(#2)).F(1))"
+          //      : "F=data(qdim(#1),qdim(#1));"
+          //      "M(#1,#2)+=sym(comp(vBase(#1).vBase(#2))(:,i,:,j).F(i,j));");
         }
         else if (normal_component) {
-          generic_assembly assem;
-          if (mf_mult.get_qdim() == 1)
-            assem.set("M(#2,#1)+=comp(Base(#2).vBase(#1).Normal())(:,:,i,i);");
-          else
-            
assem.set("M(#2,#1)+=comp(vBase(#2).mBase(#1).Normal())(:,i,:,i,j,j);");
-          assem.push_mi(mim);
-          assem.push_mf(mf_u);
-          assem.push_mf(mf_mult);
-          assem.push_mat(gmm::real_part(*B));
-          assem.assembly(rg);
+         ga_workspace workspace;
+         gmm::sub_interval Imult(0, mf_mult.nb_dof()), Iu(0, mf_u.nb_dof());
+         base_vector mult(mf_mult.nb_dof()), u(mf_u.nb_dof());
+         workspace.add_fem_variable("mult", mf_mult, Imult, mult);
+         workspace.add_fem_variable("u", mf_u, Iu, u);
+         workspace.add_expression("Test_mult.(Test2_u.Normal)", mim, rg);
+         model_real_sparse_matrix BB(mf_mult.nb_dof(), mf_u.nb_dof());
+         workspace.set_assembled_matrix(BB);
+         workspace.assembly(2);
+         gmm::add(BB, B);
+
+          // generic_assembly assem;
+          // if (mf_mult.get_qdim() == 1)
+          //   
assem.set("M(#2,#1)+=comp(Base(#2).vBase(#1).Normal())(:,:,i,i);");
+          // else
+          //   
assem.set("M(#2,#1)+=comp(vBase(#2).mBase(#1).Normal())(:,i,:,i,j,j);");
+          // assem.push_mi(mim);
+          // assem.push_mf(mf_u);
+          // assem.push_mf(mf_mult);
+          // assem.push_mat(gmm::real_part(*B));
+          // assem.assembly(rg);
         } else {
           asm_mass_matrix(*B, mim, mf_mult, mf_u, rg);
         }




reply via email to

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