[Getfem-commits] r4668 - in /trunk/getfem/interface/src: gf_asm.cc gf_compute.cc

[Yves . Renard]

Author: renard
Date: Thu May 29 11:01:10 2014
New Revision: 4668

URL: http://svn.gna.org/viewcvs/getfem?rev=4668&view=rev
Log:
allowing interpolation matrix for set of points

Modified:
    trunk/getfem/interface/src/gf_asm.cc
    trunk/getfem/interface/src/gf_compute.cc

Modified: trunk/getfem/interface/src/gf_asm.cc
URL: 
http://svn.gna.org/viewcvs/getfem/trunk/getfem/interface/src/gf_asm.cc?rev=4668&r1=4667&r2=4668&view=diff
==============================================================================
--- trunk/getfem/interface/src/gf_asm.cc        (original)
+++ trunk/getfem/interface/src/gf_asm.cc        Thu May 29 11:01:10 2014
@@ -502,12 +502,41 @@
 
 void interpolate_or_extrapolate(mexargs_in &in, mexargs_out &out, int 
extrapolate) {
   const getfem::mesh_fem *mf1 = in.pop().to_const_mesh_fem();
-  const getfem::mesh_fem *mf2 = in.pop().to_const_mesh_fem();
-  gmm::row_matrix<getfem::model_real_sparse_vector> Maux(mf2->nb_dof(), 
mf1->nb_dof());
-  getfem::interpolation(*mf1, *mf2, Maux, extrapolate);
-  gf_real_sparse_by_col M(mf2->nb_dof(), mf1->nb_dof());
-  gmm::copy(Maux, M);
-  out.pop().from_sparse(M);
+  if (in.front().is_mesh_fem()) {
+    const getfem::mesh_fem *mf2 = in.pop().to_const_mesh_fem();
+    gmm::row_matrix<getfem::model_real_sparse_vector>
+      Maux(mf2->nb_dof(), mf1->nb_dof());
+    getfem::interpolation(*mf1, *mf2, Maux, extrapolate);
+    gf_real_sparse_by_col M(mf2->nb_dof(), mf1->nb_dof());
+    gmm::copy(Maux, M);
+    out.pop().from_sparse(M);
+  } else {
+    size_type N = mf1->linked_mesh().dim();
+    darray st = in.pop().to_darray();
+    std::vector<double> PTS(st.begin(), st.end());
+    size_type nbpoints = gmm::vect_size(PTS) / N;
+    getfem::base_node p(N);
+    getfem::mesh_trans_inv mti(mf1->linked_mesh());
+    for (size_type i = 0; i < nbpoints; ++i) {
+      gmm::copy(gmm::sub_vector(PTS, gmm::sub_interval(i*N, N)), p);
+      mti.add_point(p);
+    }
+    
+    size_type qmult = mf1->get_qdim();
+    /* if (qmult != 1) dims.push_back(unsigned(qmult)); */
+/*     dims.push_back(unsigned(nbpoints)); */
+/*     dims.opt_transform_col_vect_into_row_vect(); */
+/*     garray<T> V = out.pop().create_array(dims,T()); */
+
+    gmm::row_matrix<getfem::model_real_sparse_vector>
+      Maux(nbpoints*qmult, mf1->nb_dof());
+
+    getfem::base_vector U, V;
+    getfem::interpolation(*mf1, mti, U, V, Maux, 1, extrapolate);
+    gf_real_sparse_by_col M(nbpoints*qmult, mf1->nb_dof());
+    gmm::copy(Maux, M);
+    out.pop().from_sparse(M);
+  }
 }
 
 static const getfem::mesh_im *get_mim(mexargs_in &in) {
@@ -1018,7 +1047,7 @@
       High-level generic assembly procedure for volumic assembly.
 
       Performs the generic assembly of `expression` with the integration
-      method `mim` on the mesh region of index `region' (-1 means all 
+      method `mim` on the mesh region of index `region` (-1 means all 
       the element of the mesh). The smae mesh should be shared by
       the integration method and all the finite element methods
       corresponding to the variables.
@@ -1098,8 +1127,8 @@
        );
 
 
-    /address@hidden Mi = ('interpolation matrix', @tmf mf, @tmf mfi)
-    Build the interpolation matrix from a @tmf onto another @tmf.
+    /address@hidden Mi = ('interpolation matrix', @tmf mf, address@hidden mfi 
| @vec pts} )
+    Build the interpolation matrix from a @tmf onto another @tmf or a set of 
points.
 
     Return a matrix `Mi`, such that `V = Mi.U` is equal to
     ::COMPUTE('interpolate_on',mfi). Useful for repeated interpolations.
@@ -1116,8 +1145,8 @@
        );
 
 
-    /address@hidden Me = ('extrapolation matrix',@tmf mf, @tmf mfe)
-    Build the extrapolation matrix from a @tmf onto another @tmf.
+    /address@hidden Me = ('extrapolation matrix',@tmf mf,  address@hidden mfe 
| @vec pts})
+    Build the extrapolation matrix from a @tmf onto another @tmf or a set of 
points.
 
     Return a matrix `Me`, such that `V = Me.U` is equal to
     ::COMPUTE('extrapolate_on',mfe). Useful for repeated

Modified: trunk/getfem/interface/src/gf_compute.cc
URL: 
http://svn.gna.org/viewcvs/getfem/trunk/getfem/interface/src/gf_compute.cc?rev=4668&r1=4667&r2=4668&view=diff
==============================================================================
--- trunk/getfem/interface/src/gf_compute.cc    (original)
+++ trunk/getfem/interface/src/gf_compute.cc    Thu May 29 11:01:10 2014
@@ -153,7 +153,7 @@
     getfem::mesh_trans_inv mti(mf.linked_mesh());
     for (size_type i = 0; i < nbpoints; ++i) {
       gmm::copy(gmm::sub_vector(PTS, gmm::sub_interval(i*N, N)), p);
-      cout << "adding point" << p << endl;
+      // cout << "adding point" << p << endl;
       mti.add_point(p);
     }
     
@@ -164,9 +164,9 @@
     garray<T> V = out.pop().create_array(dims,T());
 
     getfem::base_matrix Maux;
-    cout << "begin interpolation, qmult = " << qmult << endl;
+    // cout << "begin interpolation, qmult = " << qmult << endl;
     getfem::interpolation(mf, mti, U, V, Maux, 0);
-    cout << "end interpolation" << endl;
+    // cout << "end interpolation" << endl;
 
   }
   // else THROW_BADARG("expecting a mesh_fem or a mesh_slice for 
interpolation");