Merge branch 'fix_CalcElementMatrix_multidim' into 'master'

fix dimension of matrix given to CalcElementMatrix See merge request !243
NGSolve · Oct 3, 2017 · 9513c99468f3d8a24c6a0d58d65899e00833b90d · 9513c99
2 parents 3b9cbb2 + 87dddff
commit 9513c99468f3d8a24c6a0d58d65899e00833b90d
Unified Split

Showing with 33 additions and 12 deletions.

+30 −12 fem/python_fem.cpp

+3 −0 fem/symbolicintegrator.hpp
diff --git a/fem/python_fem.cpp b/fem/python_fem.cpp
@@ -1310,25 +1310,34 @@ void NGS_DLL_HEADER ExportNgfem(py::module &m) {

    .def("CalcElementMatrix",
         [] (shared_ptr<BFI> self,
-                             const FiniteElement & fe, const ElementTransformation &trafo,
-                             int heapsize)
+             const FiniteElement & fe, const ElementTransformation &trafo,
+             int heapsize, bool complex)
                         {
-                           Matrix<> mat(fe.GetNDof());
                           while (true)
                             {
                               try
                                 {
                                   LocalHeap lh(heapsize);
-                                   self->CalcElementMatrix (fe, trafo, mat, lh);
-                                   return mat;
+                                   if (complex)
+                                     {
+                                       Matrix<Complex> mat(fe.GetNDof() * self->GetDimension());
+                                       self->CalcElementMatrix(fe,trafo,mat,lh);
+                                       return py::cast(mat);
+                                     }
+                                   else
+                                     {
+                                       Matrix<> mat(fe.GetNDof() * self->GetDimension());
+                                       self->CalcElementMatrix (fe, trafo, mat, lh);
+                                       return py::cast(mat);
+                                     }
                                 }
                               catch (LocalHeapOverflow ex)
                                 {
                                   heapsize *= 10;
                                 }
-                             };
+                             }
                         },
-         py::arg("fel"),py::arg("trafo"),py::arg("heapsize")=10000)
+         py::arg("fel"),py::arg("trafo"),py::arg("heapsize")=10000, py::arg("complex") = false)
    ;


@@ -1400,24 +1409,33 @@ void NGS_DLL_HEADER ExportNgfem(py::module &m) {
        static_cast<void(LinearFormIntegrator::*)(const FiniteElement&, const ElementTransformation&, FlatVector<double>, LocalHeap&)const>(&LinearFormIntegrator::CalcElementVector))
    .def("CalcElementVector",
         [] (shared_ptr<LFI>  self, const FiniteElement & fe, const ElementTransformation& trafo,
-             int heapsize)
+             int heapsize, bool complex)
         {
-           Vector<> vec(fe.GetNDof());
           while (true)
             {
               try
                 {
                   LocalHeap lh(heapsize);
-                   self->CalcElementVector (fe, trafo, vec, lh);
-                   return vec;
+                   if (complex)
+                     {
+                       Vector<Complex> vec(fe.GetNDof() * self->GetDimension());
+                       self->CalcElementVector(fe,trafo,vec,lh);
+                       return py::cast(vec);
+                     }
+                   else
+                     {
+                       Vector<> vec(fe.GetNDof() * self->GetDimension());
+                       self->CalcElementVector (fe, trafo, vec, lh);
+                       return py::cast(vec);
+                     }
                 }
               catch (LocalHeapOverflow ex)
                 {
                   heapsize *= 10;
                 }
             };
         },
-         py::arg("fel"),py::arg("trafo"),py::arg("heapsize")=10000)
+         py::arg("fel"),py::arg("trafo"),py::arg("heapsize")=10000, py::arg("complex")=false)
    ;


diff --git a/fem/symbolicintegrator.hpp b/fem/symbolicintegrator.hpp
@@ -467,6 +467,7 @@ class CompoundDifferentialOperator : public DifferentialOperator

    virtual VorB VB() const { return vb; }
    virtual string Name () const { return string ("Symbolic LFI"); }
+    virtual int GetDimension() const override { return proxies[0]->Evaluator()->BlockDim(); }

    void SetIntegrationRule (const IntegrationRule & _ir);

@@ -523,6 +524,8 @@ class CompoundDifferentialOperator : public DifferentialOperator
    // virtual SIMD_IntegrationRule Get_SIMD_IntegrationRuleEB (const FiniteElement & fel, int facetnr, LocalHeap & lh) const;
    
    void SetIntegrationRule (const IntegrationRule & _ir);
+
+    virtual int GetDimension() const override { return trial_proxies[0]->Evaluator()->BlockDim(); }
    
    virtual void 
    CalcElementMatrix (const FiniteElement & fel,