Merge branch 'bitarray_ngcore' into 'master'

Bitarray from NGSolve to ngcore

See merge request jschoeberl/netgen!223

libsrc/core/CMakeLists.txt

@@ -1,6 +1,7 @@
 
 add_library(ngcore SHARED
   archive.cpp
+  bitarray.cpp
   localheap.cpp
   logging.cpp
   flags.cpp
@@ -58,7 +59,7 @@ target_link_libraries(ngcore PUBLIC netgen_mpi PRIVATE ${CMAKE_THREAD_LIBS_INIT}
 install(FILES ngcore.hpp archive.hpp type_traits.hpp version.hpp ngcore_api.hpp logging.hpp
   exception.hpp symboltable.hpp paje_trace.hpp utils.hpp profiler.hpp mpi_wrapper.hpp
   array.hpp taskmanager.hpp concurrentqueue.h localheap.hpp python_ngcore.hpp flags.hpp
-  xbool.hpp signal.hpp
+  xbool.hpp signal.hpp bitarray.hpp
   DESTINATION ${NG_INSTALL_DIR_INCLUDE}/core COMPONENT netgen_devel)
 
 if(ENABLE_CPP_CORE_GUIDELINES_CHECK)

libsrc/core/bitarray.cpp

@@ -0,0 +1,143 @@
+/**************************************************************************/
+/* File:   bitarray.cpp                                                   */
+/* Autho: Joachim Schoeberl                                               */
+/* Date:   01. Jun. 95                                                    */
+/**************************************************************************/
+
+/*
+   data type BitArray
+*/
+
+#include "bitarray.hpp"
+
+namespace ngcore
+{
+  BitArray :: BitArray (size_t asize)
+  {
+    size = 0;
+    data = NULL;
+    SetSize (asize);
+  }
+
+  BitArray :: BitArray (size_t asize, LocalHeap & lh)
+  {
+    size = asize;
+    data = new (lh) unsigned char [Addr (size)+1];
+    owns_data = false;
+  }
+
+  BitArray :: BitArray (const BitArray & ba2)
+  {
+    size = 0;
+    data = NULL;
+    (*this) = ba2;
+  }
+
+  void BitArray :: SetSize (size_t asize)
+  {
+    if (size == asize) return;
+    if (owns_data) delete [] data;
+
+    size = asize;
+    data = new unsigned char [Addr (size)+1];
+  }
+
+  BitArray & BitArray :: Set () throw()
+  {
+    if (!size) return *this;
+    for (size_t i = 0; i <= Addr (size); i++)
+      data[i] = UCHAR_MAX;
+    return *this;
+  }
+
+  BitArray & BitArray :: Clear () throw()
+  {
+    if (!size) return *this;
+    for (size_t i = 0; i <= Addr (size); i++)
+      data[i] = 0;
+    return *this;
+  }
+
+  BitArray & BitArray :: Invert ()
+  {
+    if (!size) return *this;
+    for (size_t i = 0; i <= Addr (size); i++)
+      data[i] ^= 255;
+    return *this;
+  }
+
+  BitArray & BitArray :: And (const BitArray & ba2)
+  {
+    if (!size) return *this;
+    for (size_t i = 0; i <= Addr (size); i++)
+      data[i] &= ba2.data[i];
+        return *this;
+  }
+
+
+  BitArray & BitArray :: Or (const BitArray & ba2)
+  {
+    if (!size) return *this;
+    for (size_t i = 0; i <= Addr (size); i++)
+      data[i] |= ba2.data[i];
+    return *this;
+  }
+
+
+  BitArray & BitArray :: operator= (const BitArray & ba2)
+  {
+    SetSize (ba2.Size());
+    if (!size)
+      return *this;
+    for (size_t i = 0; i <= Addr (size); i++)
+      data[i] = ba2.data[i];
+    return *this;
+  }
+
+  std::ostream & operator<<(std::ostream & s, const BitArray & ba)
+  {
+    size_t n = ba.Size();
+    for (size_t i = 0; i < n; i++)
+      {
+	if (i % 50 == 0) s << i << ": ";
+	s << int(ba[i]);
+	if (i % 50 == 49) s << "\n";
+      }
+    s << std::flush;
+    return s;
+  }
+
+  size_t BitArray :: NumSet () const
+  {
+    size_t cnt = 0;
+    for (size_t i = 0; i < Size(); i++)
+      if (Test(i)) cnt++;
+    return cnt;
+  }
+
+  Archive & operator & (Archive & archive, BitArray & ba)
+  {
+    if (archive.Output())
+      {
+        archive << ba.Size();
+        for (size_t i = 0; i < ba.Size(); i++)
+          archive << ba[i];
+      }
+    else
+      {
+        int size;
+        archive & size;
+        ba.SetSize (size);
+        ba.Clear();
+        for (size_t i = 0; i < size; i++)
+          {
+            bool b;
+            archive & b;
+            if (b) ba.SetBit(i);
+          }
+      }
+    return archive;
+  }
+
+
+}

libsrc/core/bitarray.hpp

@@ -0,0 +1,200 @@
+#ifndef NETGEN_CORE_BITARRAY
+#define NETGEN_CORE_BITARRAY
+
+/**************************************************************************/
+/* File:   bitarray.hpp                                                   */
+/* Author: Joachim Schoeberl                                              */
+/* Date:   01. Jun. 95                                                    */
+/**************************************************************************/
+
+#include <climits>
+#include <cstring>
+#include <ostream>
+
+#include "archive.hpp"
+#include "array.hpp"
+#include "localheap.hpp"
+#include "ngcore_api.hpp"
+#include "utils.hpp"
+
+namespace ngcore
+{
+
+/**
+   A compressed array of bools.
+
+   Provides bit-operations and whole array operations.
+*/
+class BitArray
+{
+protected:
+  /// number of bits
+  size_t size;
+
+  /// the data
+  unsigned char * data;
+  ///
+  bool owns_data = true;
+public:
+  /// empty array
+  BitArray ()
+    : size(0), data(nullptr) { ; }
+  /// array of asize bits
+  NGCORE_API BitArray (size_t asize);
+  /// array of asize bits
+  NGCORE_API BitArray (size_t asize, LocalHeap & lh);
+  ///
+  NGCORE_API BitArray (const BitArray & ba2);
+  BitArray (BitArray && ba2)
+    : size(ba2.size), data(ba2.data), owns_data(ba2.owns_data)
+  {
+    ba2.owns_data = false;
+    ba2.data = nullptr;
+  }
+
+  template <typename T>
+  NETGEN_INLINE BitArray (std::initializer_list<T> list)
+    : BitArray (list.size())
+  {
+    Clear();
+    int cnt = 0;
+    for (auto i = list.begin(); i < list.end(); i++, cnt++)
+      if (*i) SetBit(cnt);
+  }
+
+  /// delete data
+  ~BitArray ()
+  {
+    if (owns_data)
+      delete [] data;
+  }
+
+  /// Set size, loose values
+  NGCORE_API void SetSize (size_t asize);
+
+  /// the size
+  size_t Size () const { return size; }
+
+  /// set all bits
+  NGCORE_API BitArray & Set () throw();
+
+  /// clear all bits
+  NGCORE_API BitArray & Clear () throw();
+
+  /// set bit i
+  [[deprecated("Use either SetBit() or SetBitAtomic()")]]
+  void Set (size_t i) { SetBitAtomic(i); }
+
+  /// set bit i ( not thread safe )
+  void SetBit (size_t i)
+  {
+    NETGEN_CHECK_RANGE(i, 0, size);
+    data[Addr(i)] |= Mask(i);
+  }
+
+  /// set bit i ( thread safe )
+  void SetBitAtomic (size_t i)
+  {
+    NETGEN_CHECK_RANGE(i, 0, size);
+    unsigned char * p = data+Addr(i);
+    unsigned char mask = Mask(i);
+
+    AsAtomic(*p) |= mask;
+  }
+
+  /// clear bit i
+  void Clear (size_t i)
+  {
+    NETGEN_CHECK_RANGE(i, 0, size);
+    data[Addr(i)] &= ~Mask(i);
+  }
+
+  /// check bit i
+  bool Test (size_t i) const
+  {
+    NETGEN_CHECK_RANGE(i, 0, size);
+    return (data[Addr(i)] & Mask(i)) ? true : false;
+  }
+
+  /// set all bits to b
+  BitArray & operator= (bool b)
+  {
+    if (b) Set();
+    else   Clear();
+    return *this;
+  }
+
+  /// check bit i
+  bool operator[] (size_t i) const
+  {
+    NETGEN_CHECK_RANGE(i, 0, size);
+    return Test(i);
+  }
+
+
+  /// invert all bits
+  NGCORE_API BitArray & Invert ();
+
+  /// logical AND with ba2
+  NGCORE_API BitArray & And (const BitArray & ba2);
+
+  /// logical OR with ba2
+  NGCORE_API BitArray & Or (const BitArray & ba2);
+
+  /// copy from ba2
+  NGCORE_API BitArray & operator= (const BitArray & ba2);
+
+  NGCORE_API size_t NumSet () const;
+private:
+  ///
+  unsigned char Mask (size_t i) const
+  { return char(1) << (i % CHAR_BIT); }
+
+  ///
+  size_t Addr (size_t i) const
+  { return (i / CHAR_BIT); }
+
+};
+
+
+  inline BitArray & operator|= (BitArray & me, const BitArray & you)
+  {
+    me.Or(you);
+    return me;
+  }
+
+  inline BitArray & operator&= (BitArray & me, const BitArray & you)
+  {
+    me.And(you);
+    return me;
+  }
+
+  inline BitArray operator| (const BitArray & a, const BitArray & b)
+  {
+    BitArray res = a;
+    res |= b;
+    return res;
+  }
+
+  inline BitArray operator& (const BitArray & a, const BitArray & b)
+  {
+    BitArray res = a;
+    res &= b;
+    return res;
+  }
+
+  inline BitArray operator~ (const BitArray & a)
+  {
+    BitArray res = a;
+    res.Invert();
+    return res;
+  }
+
+
+  NGCORE_API std::ostream & operator<<(std::ostream & s, const BitArray & ba);
+
+  NGCORE_API Archive & operator & (Archive & archive, BitArray & ba);
+
+} // namespace ngcore
+
+#endif // NETGEN_CORE_BITARRAY

libsrc/core/ngcore.hpp

@@ -3,6 +3,7 @@
 
 #include "archive.hpp"
 #include "array.hpp"
+#include "bitarray.hpp"
 #include "exception.hpp"
 #include "flags.hpp"
 #include "localheap.hpp"

libsrc/core/python_ngcore.hpp

@@ -3,6 +3,7 @@
 
 #include "ngcore_api.hpp" // for operator new
 #include <pybind11/pybind11.h>
+#include <pybind11/operators.h>
 
 #include "array.hpp"
 #include "archive.hpp"

libsrc/core/python_ngcore_export.cpp

@@ -1,5 +1,6 @@
-
 #include "python_ngcore.hpp"
+#include "bitarray.hpp"
+
 
 using namespace ngcore;
 using namespace std;
@@ -11,6 +12,83 @@ PYBIND11_MODULE(pyngcore, m) // NOLINT
   ExportArray<size_t>(m);
   ExportArray<double>(m);
 
+  py::class_<BitArray, shared_ptr<BitArray>> (m, "BitArray")
+    .def(py::init([] (size_t n) { return make_shared<BitArray>(n); }),py::arg("n"))
+    .def(py::init([] (const BitArray& a) { return make_shared<BitArray>(a); } ), py::arg("ba"))
+    .def(py::init([] (const vector<bool> & a)
+                  {
+                    auto ba = make_shared<BitArray>(a.size());
+                    ba->Clear();
+                    for (size_t i = 0; i < a.size(); i++)
+                      if (a[i]) ba->SetBit(i);
+                    return ba;
+                  } ), py::arg("vec"))
+    .def("__str__", &ToString<BitArray>)
+    .def("__len__", &BitArray::Size)
+    .def("__getitem__", [] (BitArray & self, int i)
+                                         {
+                                           if (i < 0 || i >= self.Size())
+                                             throw py::index_error();
+                                           return self.Test(i);
+                                         }, py::arg("pos"), "Returns bit from given position")
+    .def("__setitem__", [] (BitArray & self, int i, bool b)
+                                         {
+                                           if (i < 0 || i >= self.Size())
+                                             throw py::index_error();
+                                           if (b) self.SetBit(i); else self.Clear(i);
+                                         }, py::arg("pos"), py::arg("value"), "Clear/Set bit at given position")
+
+    .def("__setitem__", [] (BitArray & self, py::slice inds, bool b)
+                                         {
+                                           size_t start, step, stop, n;
+                                           if (!inds.compute(self.Size(), &start, &stop, &step, &n))
+                                             throw py::error_already_set();
+
+                                           if (start == 0 && n == self.Size() && step == 1)
+                                             { // base branch
+                                               if (b)
+                                                 self.Set();
+                                               else
+                                                 self.Clear();
+                                             }
+                                           else
+                                             {
+                                               if (b)
+                                                 for (size_t i=0; i<n; i++, start+=step)
+                                                   self.SetBit(start);
+                                               else
+                                                 for (size_t i=0; i<n; i++, start+=step)
+                                                   self.Clear(start);
+                                             }
+                                         }, py::arg("inds"), py::arg("value"), "Clear/Set bit at given positions")
+
+    .def("__setitem__", [](BitArray & self,  IntRange range, bool b)
+      {
+        if (b)
+          for (size_t i : range)
+            self.SetBit(i);
+        else
+          for (size_t i : range)
+            self.Clear(i);
+      }, py::arg("range"), py::arg("value"), "Set value for range of indices" )
+
+    .def("NumSet", &BitArray::NumSet)
+    .def("Set", [] (BitArray & self) { self.Set(); }, "Set all bits")
+    .def("Set", &BitArray::SetBit, py::arg("i"), "Set bit at given position")
+    .def("Clear", [] (BitArray & self) { self.Clear(); }, "Clear all bits")
+    .def("Clear", [] (BitArray & self, int i)
+                                   {
+                                       self.Clear(i);
+                                   }, py::arg("i"), "Clear bit at given position")
+
+
+    .def(py::self | py::self)
+    .def(py::self & py::self)
+    .def(py::self |= py::self)
+    .def(py::self &= py::self)
+    .def(~py::self)
+    ;
+
   py::class_<Flags>(m, "Flags")
     .def(py::init<>())
     .def("__str__", &ToString<Flags>)

libsrc/csg/edgeflw.cpp

@@ -1233,7 +1233,7 @@ namespace netgen
 	*testout << "inv: " << endl << refedgesinv << endl;
       }
     
-    BitArray todelete(refedges.Size());
+    NgBitArray todelete(refedges.Size());
     todelete.Clear();
 
 
@@ -1748,7 +1748,7 @@ namespace netgen
     int nsol = geometry.GetNTopLevelObjects();
     
 
-    BitArray pointatsurface (nsurf);
+    NgBitArray pointatsurface (nsurf);
     pointatsurface.Clear();
   
     for (int i = 1; i <= mesh.GetNSeg(); i++)

libsrc/csg/specpoin.cpp

@@ -2060,7 +2060,7 @@ namespace netgen
       }
 
     /*
-    BitArray testuncond (specpoints.Size());
+    NgBitArray testuncond (specpoints.Size());
     testuncond.Clear();
     for(int i = 0; i<specpoints.Size(); i++)
       {
@@ -2093,7 +2093,7 @@ namespace netgen
     // if special point is unconditional on some solid,
     // it must be unconditional everywhere:
 
-    BitArray uncond (apoints.Size());
+    NgBitArray uncond (apoints.Size());
     uncond.Clear();
 
     for (int i = 0; i < specpoints.Size(); i++)

libsrc/csg/triapprox.cpp

@@ -28,7 +28,7 @@ namespace netgen
 
   void TriangleApproximation :: RemoveUnusedPoints ()
   {
-    BitArray used(GetNP());
+    NgBitArray used(GetNP());
     NgArray<int> map (GetNP());
     int i, j;
     int cnt = 0;

libsrc/csg/zrefine.cpp

@@ -255,7 +255,7 @@ namespace netgen
     NgArray<INDEX_3> ref_singular;
     NgArray<INDEX_4 > ref_slices;
 
-    BitArray first_id(geom->identifications.Size());
+    NgBitArray first_id(geom->identifications.Size());
     first_id.Set();
 
   

libsrc/general/CMakeLists.txt

@@ -2,14 +2,14 @@ add_definitions(-DNGINTERFACE_EXPORTS)
 add_library(gen INTERFACE)
 set(sdir ${CMAKE_CURRENT_SOURCE_DIR})
 target_sources(gen INTERFACE
-  ${sdir}/ngarray.cpp ${sdir}/bitarray.cpp ${sdir}/dynamicmem.cpp
+  ${sdir}/ngarray.cpp ${sdir}/ngbitarray.cpp ${sdir}/dynamicmem.cpp
         ${sdir}/hashtabl.cpp ${sdir}/mystring.cpp ${sdir}/optmem.cpp ${sdir}/parthreads.cpp
         ${sdir}/seti.cpp ${sdir}/sort.cpp ${sdir}/spbita2d.cpp ${sdir}/table.cpp
         ${sdir}/mpi_interface.cpp ${sdir}/gzstream.cpp
         )
 
 install(FILES
-  ngarray.hpp autodiff.hpp autoptr.hpp bitarray.hpp
+  ngarray.hpp autodiff.hpp autoptr.hpp ngbitarray.hpp
   dynamicmem.hpp hashtabl.hpp mpi_interface.hpp myadt.hpp
   ngsimd.hpp mystring.hpp netgenout.hpp ngpython.hpp
   optmem.hpp parthreads.hpp seti.hpp sort.hpp

libsrc/general/bitarray.hpp

@@ -1,227 +0,0 @@
-#ifndef FILE_BitArray
-#define FILE_BitArray
-
-/**************************************************************************/
-/* File:   bitarray.hpp                                                   */
-/* Author: Joachim Schoeberl                                              */
-/* Date:   01. Jun. 95                                                    */
-/**************************************************************************/
-
-#include <limits.h>
-
-namespace netgen
-{
-
-
-/**
-   data type BitArray
-   
-   BitArray is a compressed array of Boolean information. By Set and Clear
-   the whole array or one bit can be set or reset, respectively. 
-   Test returns the state of the occurring bit.
-   No range checking is done.
-
-   index ranges from 0 to size-1
-*/
-class BitArray
-{
-  INDEX size;
-  unsigned char * data;
-
-public:
-  BitArray ();
-  ///
-  BitArray (INDEX asize);
-  ///
-  ~BitArray ();
-
-  /// 
-  void SetSize (INDEX asize);
-  ///
-  INDEX Size () const
-  {
-    return size;
-  }
-
-  ///
-  void Set ();
-  ///
-  void Set (INDEX i)
-  {
-    data[Addr(i)] |= Mask(i);
-  }
-  
-  void Clear ();
-
-
-  void Clear (INDEX i)
-  {
-    data[Addr(i)] &= ~Mask(i);
-  }
-
-  bool Test (INDEX i) const
-  {
-    return (data[i / CHAR_BIT] & (char(1) << (i % CHAR_BIT) ) ) ? true : false;
-  }
-
-  ///
-  void Invert ();
-  ///
-  void And (const BitArray & ba2);
-  ///
-  void Or (const BitArray & ba2);
-private:
-  ///
-  inline unsigned char Mask (INDEX i) const
-  {
-    return char(1) << (i % CHAR_BIT);
-  }
-  ///
-  inline INDEX Addr (INDEX i) const
-  {
-  return (i / CHAR_BIT);
-  }
-
-  ///
-  BitArray & operator= (BitArray &);
-  ///
-  BitArray (const BitArray &);
-};
-
-
-
-// print bitarray
-inline ostream & operator<< (ostream & s, const BitArray & a)
-{
-  for (int i = 1; i <= a.Size(); i++)
-    {
-      s << int (a.Test(i));
-      if (i % 40 == 0) s << "\n";
-    }
-  if (a.Size() % 40 != 0) s << "\n";
-  return s;
-}
-
-
-/*
-inline
-INDEX BitArray :: Size () const
-  {
-  return size;
-  }
-
-inline
-unsigned char BitArray :: Mask (INDEX i) const
-  {
-  return char(1) << (i % CHAR_BIT);
-  }
-
-inline
-INDEX BitArray :: Addr (INDEX i) const
-  {
-  return (i / CHAR_BIT);
-  }
-inline
-void BitArray :: Set (INDEX i)
-  {
-  data[Addr(i)] |= Mask(i);
-  }
-
-inline
-void BitArray :: Clear (INDEX i)
-  {
-  data[Addr(i)] &= ~Mask(i);
-  }
-
-
-inline
-int BitArray :: Test (INDEX i) const
-  {
-  return (data[i / CHAR_BIT] & (char(1) << (i % CHAR_BIT) ) ) ? 1 : 0;
-  }
-
-*/
-
-
-
-
-
-
-/**
-   data type BitArrayChar
-   
-   BitArray is an array of Boolean information. By Set and Clear
-   the whole array or one bit can be set or reset, respectively. 
-   Test returns the state of the occurring bit.
-   No range checking is done.
-*/
-template <int BASE = 1>
-class BitArrayChar
-{
-  ///
-  NgArray<char,BASE> data;
-
-public:
-  ///
-  BitArrayChar ()
-  { ; }
-  ///
-  BitArrayChar (int asize)
-    : data(asize)
-  { ; }
-  ///
-  ~BitArrayChar ()
-  { ; }
-
-  ///
-  void SetSize (int asize)
-  { data.SetSize(asize); }
-
-  ///
-  inline int Size () const
-  { return data.Size(); }
-
-  ///
-  void Set ();
-  ///
-  inline void Set (int i)
-  { data[i] = 1; }
-  ///
-  void Clear ();
-  ///
-  inline void Clear (int i)
-  { data[i] = 0; }
-  ///
-  inline int Test (int i) const
-  { return data[i]; }
-  ///
-  void Invert ();
-  ///
-  void And (const BitArrayChar & ba2);
-  ///
-  void Or (const BitArrayChar & ba2);
-private:
-  ///  copy bitarray is not supported
-  BitArrayChar & operator= (BitArrayChar &) { return *this; }
-  ///  copy bitarray is not supported
-  BitArrayChar (const BitArrayChar &) { ; }
-};
-
-
-
-
-template <int BASE>
-inline ostream & operator<< (ostream & s, const BitArrayChar<BASE> & a)
-{
-  for (int i = BASE; i < a.Size()+BASE; i++)
-    {
-      s << a.Test(i);
-      if ( (i-BASE) % 40 == 39) s << "\n";
-    }
-  if (a.Size() % 40 != 0) s << "\n";
-  return s;
-}
-
-}
-
-#endif

libsrc/general/myadt.hpp

@@ -33,7 +33,7 @@ namespace netgen
 #include "hashtabl.hpp"
 
 
-#include "bitarray.hpp"
+#include "ngbitarray.hpp"
 #include "spbita2d.hpp"
 
 #include "seti.hpp"

libsrc/general/ngbitarray.cpp

@@ -5,7 +5,7 @@
 /**************************************************************************/
 
 /* 
-   data type BitArray
+   data type NgBitArray
 */
 
 #include <mystdlib.h>
@@ -16,25 +16,25 @@ namespace netgen
 {
   //using namespace netgen;
 
-  BitArray :: BitArray ()
+  NgBitArray :: NgBitArray ()
   {
     size = 0;
     data = NULL;
   }
 
-  BitArray :: BitArray (int asize)
+  NgBitArray :: NgBitArray (int asize)
   {
     size = 0;
     data = NULL;
     SetSize (asize);
   }
 
-  BitArray :: ~BitArray ()
+  NgBitArray :: ~NgBitArray ()
   {
     delete [] data;
   }
 
-  void BitArray :: SetSize (int asize)
+  void NgBitArray :: SetSize (int asize)
   {
     if (size == asize) return;
     delete [] data;
@@ -43,14 +43,14 @@ namespace netgen
     data = new unsigned char [Addr (size)+1];
   }
 
-  void BitArray :: Set ()
+  void NgBitArray :: Set ()
   {
     if (!size) return;
     for (int i = 0; i <= Addr (size); i++)
       data[i] = UCHAR_MAX;
   }
 
-  void BitArray :: Clear ()
+  void NgBitArray :: Clear ()
   {
     if (!size) return;
     for (int i = 0; i <= Addr (size); i++)
@@ -59,14 +59,14 @@ namespace netgen
 
 
 
-  void BitArray :: Invert ()
+  void NgBitArray :: Invert ()
   {
     if (!size) return;
     for (int i = 0; i <= Addr (size); i++)
       data[i] ^= 255;
   }
 
-  void BitArray :: And (const BitArray & ba2)
+  void NgBitArray :: And (const NgBitArray & ba2)
   {
     if (!size) return;
     for (int i = 0; i <= Addr (size); i++)
@@ -74,7 +74,7 @@ namespace netgen
   }
 
 
-  void BitArray :: Or (const BitArray & ba2)
+  void NgBitArray :: Or (const NgBitArray & ba2)
   {
     if (!size) return;
     for (int i = 0; i <= Addr (size); i++)
@@ -82,51 +82,4 @@ namespace netgen
   }
 
 
-
-
-
-
-
-
-
-
-
-  template <int BASE>
-  void BitArrayChar<BASE> :: Set ()
-  {
-    data = 1;
-  }
-
-  template <int BASE>
-  void BitArrayChar<BASE> :: Clear ()
-  {
-    data = 0;
-  }
-
-
-  template <int BASE>
-  void BitArrayChar<BASE> :: Invert ()
-  {
-    for (int i = BASE; i < data.Size()+BASE; i++)
-      data[i] = 1 - data[i];
-  }
-
-  template <int BASE>
-  void BitArrayChar<BASE> :: And (const BitArrayChar & ba2)
-  {
-    for (int i = BASE; i < data.Size()+BASE; i++)
-      data[i] &= ba2.data[i];
-  }
-  
-
-  template <int BASE>
-  void BitArrayChar<BASE> :: Or (const BitArrayChar & ba2)
-  {
-    for (int i = BASE; i < data.Size()+BASE; i++)
-      data[i] |= ba2.data[i];
-  }
-  
-
-  template class BitArrayChar<0>;
-  template class BitArrayChar<1>;
 }

libsrc/general/ngbitarray.hpp

@@ -0,0 +1,147 @@
+#ifndef FILE_BitArray
+#define FILE_BitArray
+
+/**************************************************************************/
+/* File:   bitarray.hpp                                                   */
+/* Author: Joachim Schoeberl                                              */
+/* Date:   01. Jun. 95                                                    */
+/**************************************************************************/
+
+#include <limits.h>
+
+namespace netgen
+{
+
+
+/**
+   data type NgBitArray
+   
+   NgBitArray is a compressed array of Boolean information. By Set and Clear
+   the whole array or one bit can be set or reset, respectively. 
+   Test returns the state of the occurring bit.
+   No range checking is done.
+
+   index ranges from 0 to size-1
+*/
+class NgBitArray
+{
+  INDEX size;
+  unsigned char * data;
+
+public:
+  NgBitArray ();
+  ///
+  NgBitArray (INDEX asize);
+  ///
+  ~NgBitArray ();
+
+  /// 
+  void SetSize (INDEX asize);
+  ///
+  INDEX Size () const
+  {
+    return size;
+  }
+
+  ///
+  void Set ();
+  ///
+  void Set (INDEX i)
+  {
+    data[Addr(i)] |= Mask(i);
+  }
+  
+  void Clear ();
+
+
+  void Clear (INDEX i)
+  {
+    data[Addr(i)] &= ~Mask(i);
+  }
+
+  bool Test (INDEX i) const
+  {
+    return (data[i / CHAR_BIT] & (char(1) << (i % CHAR_BIT) ) ) ? true : false;
+  }
+
+  ///
+  void Invert ();
+  ///
+  void And (const NgBitArray & ba2);
+  ///
+  void Or (const NgBitArray & ba2);
+private:
+  ///
+  inline unsigned char Mask (INDEX i) const
+  {
+    return char(1) << (i % CHAR_BIT);
+  }
+  ///
+  inline INDEX Addr (INDEX i) const
+  {
+  return (i / CHAR_BIT);
+  }
+
+  ///
+  NgBitArray & operator= (NgBitArray &);
+  ///
+  NgBitArray (const NgBitArray &);
+};
+
+
+
+// print bitarray
+inline ostream & operator<< (ostream & s, const NgBitArray & a)
+{
+  for (int i = 1; i <= a.Size(); i++)
+    {
+      s << int (a.Test(i));
+      if (i % 40 == 0) s << "\n";
+    }
+  if (a.Size() % 40 != 0) s << "\n";
+  return s;
+}
+
+
+/*
+inline
+INDEX NgBitArray :: Size () const
+  {
+  return size;
+  }
+
+inline
+unsigned char NgBitArray :: Mask (INDEX i) const
+  {
+  return char(1) << (i % CHAR_BIT);
+  }
+
+inline
+INDEX NgBitArray :: Addr (INDEX i) const
+  {
+  return (i / CHAR_BIT);
+  }
+inline
+void NgBitArray :: Set (INDEX i)
+  {
+  data[Addr(i)] |= Mask(i);
+  }
+
+inline
+void NgBitArray :: Clear (INDEX i)
+  {
+  data[Addr(i)] &= ~Mask(i);
+  }
+
+
+inline
+int NgBitArray :: Test (INDEX i) const
+  {
+  return (data[i / CHAR_BIT] & (char(1) << (i % CHAR_BIT) ) ) ? 1 : 0;
+  }
+
+*/
+
+}
+
+#endif

libsrc/general/seti.hpp

@@ -17,7 +17,7 @@ namespace netgen
 class IndexSet
 {
   NgArray<int> set;
-  BitArray flags;
+  NgBitArray flags;
 public:
   IndexSet (int maxind);
   

libsrc/interface/writeabaqus.cpp

@@ -140,7 +140,7 @@ void WriteAbaqusFormat (const Mesh & mesh,
       int masternode(0);
 
       NgArray<INDEX_2> pairs;
-      BitArray master(np), help(np);
+      NgBitArray master(np), help(np);
       master.Set();
       for (i = 1; i <= 3; i++)
 	{
@@ -205,7 +205,7 @@ void WriteAbaqusFormat (const Mesh & mesh,
 	      << "*EQUATION, INPUT=" << mpcfilename << endl;
 	  
 
-      BitArray eliminated(np);
+      NgBitArray eliminated(np);
       eliminated.Clear();
       for (i = 1; i <= mesh.GetIdentifications().GetMaxNr(); i++)
 	{

libsrc/interface/writeuser.cpp

@@ -938,7 +938,7 @@ void WriteFile (int typ,
       NgArray<INDEX_2> edgelist;
 
       // edge (point) on boundary ?
-      BitArray bedge, bpoint(mesh.GetNP());
+      NgBitArray bedge, bpoint(mesh.GetNP());
 
       static int eledges[6][2] = { { 1, 2 } , { 1, 3 } , { 1, 4 },
 				   { 2, 3 } , { 2, 4 } , { 3, 4 } };

libsrc/meshing/adfront3.cpp

@@ -349,18 +349,18 @@ void AdFront3 :: RebuildInternalTables ()
 
 
 
-  BitArrayChar<PointIndex::BASE> usecl(np);
+  Array<bool, PointIndex> usecl(np);
-  usecl.Clear();
+  usecl = false;
   for (int i = 1; i <= faces.Size(); i++)
     {
-      usecl.Set (points[faces.Get(i).Face().PNum(1)].cluster);
+      usecl[points[faces.Get(i).Face().PNum(1)].cluster] = true;
       faces.Elem(i).cluster =
 	points[faces.Get(i).Face().PNum(1)].cluster;
     }
   int cntcl = 0;
   for (int i = PointIndex::BASE; 
        i < np+PointIndex::BASE; i++)
-    if (usecl.Test(i))
+    if (usecl[i])
       cntcl++;
 
   NgArray<double, PointIndex::BASE> clvol (np);

libsrc/meshing/bisect.cpp

@@ -3159,7 +3159,7 @@ namespace netgen
 	if (opt.refine_hp)
 	  {
 	    PrintMessage(3,"refine hp");
-	    BitArray singv(np);
+	    NgBitArray singv(np);
 	    singv.Clear();
 
 	    if (mesh.GetDimension() == 3)
@@ -3833,7 +3833,7 @@ namespace netgen
       }
     */
 
-    BitArray isnewpoint(np);
+    NgBitArray isnewpoint(np);
     isnewpoint.Clear();
 
     for (int i = 0; i < cutedges.Size(); i++)

libsrc/meshing/boundarylayer.cpp

@@ -19,7 +19,7 @@ namespace netgen
       cout << "Old NP: " << mesh.GetNP() << endl;
       cout << "Trigs: " << mesh.GetNSE() << endl;
 
-      BitArray bndnodes(np);
+      NgBitArray bndnodes(np);
       NgArray<int> mapto(np);
 
       bndnodes.Clear();
@@ -210,7 +210,7 @@ namespace netgen
          int nseg = mesh.GetNSeg();
 
          // Indicate which points need to be remapped
-         BitArray bndnodes(np+1);  // big enough for 1-based array
+         NgBitArray bndnodes(np+1);  // big enough for 1-based array
 
          // Map of the old points to the new points
          NgArray<PointIndex, PointIndex::BASE> mapto(np);
@@ -286,7 +286,7 @@ namespace netgen
          // don't have boundary layers
 
          // Bit array to keep track of segments already processed
-         BitArray segsel(nseg);
+         NgBitArray segsel(nseg);
 
          // Set them all to "1" to initially activate all segments
          segsel.Set();

libsrc/meshing/classifyhpel.hpp

@@ -1,5 +1,5 @@
 HPREF_ELEMENT_TYPE ClassifyTet(HPRefElement & el, INDEX_2_HASHTABLE<int> & edges, INDEX_2_HASHTABLE<int> & edgepoint_dom, 
-                               BitArray & cornerpoint, BitArray & edgepoint, INDEX_3_HASHTABLE<int> & faces, INDEX_2_HASHTABLE<int> & face_edges, 
+                               NgBitArray & cornerpoint, NgBitArray & edgepoint, INDEX_3_HASHTABLE<int> & faces, INDEX_2_HASHTABLE<int> & face_edges, 
                                INDEX_2_HASHTABLE<int> & surf_edges, NgArray<int, PointIndex::BASE> & facepoint)
 {
   int ep1(0), ep2(0), ep3(0), ep4(0), cp1(0), cp2(0), cp3(0), cp4(0), fp1, fp2, fp3, fp4;
@@ -423,7 +423,7 @@ HPREF_ELEMENT_TYPE ClassifyTet(HPRefElement & el, INDEX_2_HASHTABLE<int> & edges
 
 
 HPREF_ELEMENT_TYPE ClassifyPrism(HPRefElement & el, INDEX_2_HASHTABLE<int> & edges, INDEX_2_HASHTABLE<int> & edgepoint_dom, 
-                                 BitArray & cornerpoint, BitArray & edgepoint, INDEX_3_HASHTABLE<int> & faces, INDEX_2_HASHTABLE<int> & face_edges, 
+                                 NgBitArray & cornerpoint, NgBitArray & edgepoint, INDEX_3_HASHTABLE<int> & faces, INDEX_2_HASHTABLE<int> & face_edges, 
                                  INDEX_2_HASHTABLE<int> & surf_edges, NgArray<int, PointIndex::BASE> & facepoint)
 {
 
@@ -660,7 +660,7 @@ HPREF_ELEMENT_TYPE ClassifyPrism(HPRefElement & el, INDEX_2_HASHTABLE<int> & edg
 
 // #ifdef SABINE 
 HPREF_ELEMENT_TYPE ClassifyTrig(HPRefElement & el, INDEX_2_HASHTABLE<int> & edges, INDEX_2_HASHTABLE<int> & edgepoint_dom, 
-                                BitArray & cornerpoint, BitArray & edgepoint, INDEX_3_HASHTABLE<int> & faces, INDEX_2_HASHTABLE<int>  & face_edges, 
+                                NgBitArray & cornerpoint, NgBitArray & edgepoint, INDEX_3_HASHTABLE<int> & faces, INDEX_2_HASHTABLE<int>  & face_edges, 
 				INDEX_2_HASHTABLE<int> & surf_edges, NgArray<int, PointIndex::BASE> & facepoint, int dim, const FaceDescriptor & fd)
 
 {
@@ -875,7 +875,7 @@ HPREF_ELEMENT_TYPE ClassifyTrig(HPRefElement & el, INDEX_2_HASHTABLE<int> & edge
 }
 #ifdef HPREF_OLD 
 HPREF_ELEMENT_TYPE ClassifyTrig(HPRefElement & el, INDEX_2_HASHTABLE<int> & edges, INDEX_2_HASHTABLE<int> & edgepoint_dom, 
-				BitArray & cornerpoint, BitArray & edgepoint, INDEX_3_HASHTABLE<int> & faces, INDEX_2_HASHTABLE<int> & face_edges, 
+				NgBitArray & cornerpoint, NgBitArray & edgepoint, INDEX_3_HASHTABLE<int> & faces, INDEX_2_HASHTABLE<int> & face_edges, 
 				INDEX_2_HASHTABLE<int> & surf_edges, NgArray<int, PointIndex::BASE> & facepoint, int dim, const FaceDescriptor & fd)
 {
   HPREF_ELEMENT_TYPE type = HP_NONE; 
@@ -1136,7 +1136,7 @@ HPREF_ELEMENT_TYPE ClassifyTrig(HPRefElement & el, INDEX_2_HASHTABLE<int> & edge
 }
 #endif
 HPREF_ELEMENT_TYPE ClassifyQuad(HPRefElement & el, INDEX_2_HASHTABLE<int> & edges, INDEX_2_HASHTABLE<int> & edgepoint_dom, 
-                                BitArray & cornerpoint, BitArray & edgepoint, INDEX_3_HASHTABLE<int> & faces, INDEX_2_HASHTABLE<int>  & face_edges, 
+                                NgBitArray & cornerpoint, NgBitArray & edgepoint, INDEX_3_HASHTABLE<int> & faces, INDEX_2_HASHTABLE<int>  & face_edges, 
 				INDEX_2_HASHTABLE<int> & surf_edges, NgArray<int, PointIndex::BASE> & facepoint, int dim, const FaceDescriptor & fd)
 {
   HPREF_ELEMENT_TYPE type = HP_NONE; 
@@ -1486,7 +1486,7 @@ HPREF_ELEMENT_TYPE ClassifyQuad(HPRefElement & el, INDEX_2_HASHTABLE<int> & edge
 
 
 HPREF_ELEMENT_TYPE ClassifyHex(HPRefElement & el, INDEX_2_HASHTABLE<int> & edges, INDEX_2_HASHTABLE<int> & edgepoint_dom, 
-                               BitArray & cornerpoint, BitArray & edgepoint, INDEX_3_HASHTABLE<int> & faces, INDEX_2_HASHTABLE<int> & face_edges, 
+                               NgBitArray & cornerpoint, NgBitArray & edgepoint, INDEX_3_HASHTABLE<int> & faces, INDEX_2_HASHTABLE<int> & face_edges, 
                                INDEX_2_HASHTABLE<int> & surf_edges, NgArray<int, PointIndex::BASE> & facepoint)
 {
   HPREF_ELEMENT_TYPE type = HP_NONE;
@@ -1586,7 +1586,7 @@ HPREF_ELEMENT_TYPE ClassifyHex(HPRefElement & el, INDEX_2_HASHTABLE<int> & edges
 }
 
 HPREF_ELEMENT_TYPE ClassifySegm(HPRefElement & hpel, INDEX_2_HASHTABLE<int> & edges, INDEX_2_HASHTABLE<int> & edgepoint_dom, 
-                                BitArray & cornerpoint, BitArray & edgepoint, INDEX_3_HASHTABLE<int> & faces, INDEX_2_HASHTABLE<int> & face_edges, 
+                                NgBitArray & cornerpoint, NgBitArray & edgepoint, INDEX_3_HASHTABLE<int> & faces, INDEX_2_HASHTABLE<int> & face_edges, 
                                 INDEX_2_HASHTABLE<int> & surf_edges, NgArray<int, PointIndex::BASE> & facepoint)
 {
   
@@ -1629,7 +1629,7 @@ HPREF_ELEMENT_TYPE ClassifySegm(HPRefElement & hpel, INDEX_2_HASHTABLE<int> & ed
 
 
 HPREF_ELEMENT_TYPE ClassifyPyramid(HPRefElement & el, INDEX_2_HASHTABLE<int> & edges, INDEX_2_HASHTABLE<int> & edgepoint_dom, 
-                                   BitArray & cornerpoint, BitArray & edgepoint, INDEX_3_HASHTABLE<int> & faces, INDEX_2_HASHTABLE<int> & face_edges, 
+                                   NgBitArray & cornerpoint, NgBitArray & edgepoint, INDEX_3_HASHTABLE<int> & faces, INDEX_2_HASHTABLE<int> & face_edges, 
                                    INDEX_2_HASHTABLE<int> & surf_edges, NgArray<int, PointIndex::BASE> & facepoint)
 {
   HPREF_ELEMENT_TYPE type = HP_NONE;

libsrc/meshing/delaunay.cpp

@@ -561,19 +561,19 @@ namespace netgen
     startel[3] = mesh.AddPoint (cp4);
 
     // flag points to use for Delaunay:
-    BitArrayChar<PointIndex::BASE> usep(np);
+    Array<bool, PointIndex> usep(np);
-    usep.Clear();
+    usep = false;
 
     for (auto & face : adfront->Faces())
       for (PointIndex pi : face.Face().PNums())      
-        usep.Set (pi);
+        usep[pi] = true;
     
     for (size_t i = oldnp + PointIndex::BASE; 
 	 i < np + PointIndex::BASE; i++)
-      usep.Set (i);
+      usep[i] = true;
 
     for (PointIndex pi : mesh.LockedPoints())
-      usep.Set (pi);
+      usep[pi] = true;
     
 
     NgArray<int> freelist;
@@ -649,7 +649,7 @@ namespace netgen
 
 	PointIndex newpi = mixed[pi];
 
-	if (!usep.Test(newpi)) 
+	if (!usep[newpi])
 	  continue;
 
 	cntp++;
@@ -805,7 +805,7 @@ namespace netgen
 
     // remove degenerated
 
-    BitArray badnode(mesh.GetNP());
+    NgBitArray badnode(mesh.GetNP());
     badnode.Clear();
     int ndeg = 0;
     for (int i = 1; i <= tempels.Size(); i++)
@@ -1320,7 +1320,7 @@ namespace netgen
  
 
     ne = tempels.Size();
-    BitArray inner(ne), outer(ne);
+    NgBitArray inner(ne), outer(ne);
     inner.Clear();
     outer.Clear();
     NgArray<int> elstack;

libsrc/meshing/hprefinement.cpp

@@ -551,7 +551,7 @@ namespace netgen
   }
 
   bool CheckSingularities(Mesh & mesh, INDEX_2_HASHTABLE<int> & edges, INDEX_2_HASHTABLE<int> & edgepoiclt_dom, 
-		       BitArray & cornerpoint, BitArray & edgepoint, INDEX_3_HASHTABLE<int> & faces, INDEX_2_HASHTABLE<int> & face_edges, 
+		       NgBitArray & cornerpoint, NgBitArray & edgepoint, INDEX_3_HASHTABLE<int> & faces, INDEX_2_HASHTABLE<int> & face_edges, 
 			INDEX_2_HASHTABLE<int> & surf_edges, NgArray<int, PointIndex::BASE> & facepoint, int & levels, int & act_ref); 
 
   bool ClassifyHPElements (Mesh & mesh, NgArray<HPRefElement> & elements, int & act_ref, int & levels);
@@ -1558,7 +1558,7 @@ namespace netgen
   }
 
 bool CheckSingularities(Mesh & mesh, INDEX_2_HASHTABLE<int> & edges, INDEX_2_HASHTABLE<int> & edgepoint_dom, 
-		       BitArray & cornerpoint, BitArray & edgepoint, INDEX_3_HASHTABLE<int> & faces, INDEX_2_HASHTABLE<int> & face_edges, 
+		       NgBitArray & cornerpoint, NgBitArray & edgepoint, INDEX_3_HASHTABLE<int> & faces, INDEX_2_HASHTABLE<int> & face_edges, 
 			INDEX_2_HASHTABLE<int> & surf_edges, NgArray<int, PointIndex::BASE> & facepoint, int & levels, int & act_ref)
 { 
   bool sing = 0; 
@@ -1812,11 +1812,11 @@ bool CheckSingularities(Mesh & mesh, INDEX_2_HASHTABLE<int> & edges, INDEX_2_HAS
   bool ClassifyHPElements (Mesh & mesh, NgArray<HPRefElement> & elements, int & act_ref, int & levels)
   {
     INDEX_2_HASHTABLE<int> edges(mesh.GetNSeg()+1);
-    BitArray edgepoint(mesh.GetNP());
+    NgBitArray edgepoint(mesh.GetNP());
     INDEX_2_HASHTABLE<int> edgepoint_dom(mesh.GetNSeg()+1);
 
     edgepoint.Clear();
-    BitArray cornerpoint(mesh.GetNP());
+    NgBitArray cornerpoint(mesh.GetNP());
     cornerpoint.Clear();
 
     // value = nr > 0 ... refine elements in domain nr

libsrc/meshing/improve3.cpp

@@ -607,11 +607,11 @@ void MeshOptimize3d :: SplitImprove (Mesh & mesh,
   TABLE<ElementIndex,PointIndex::BASE> elementsonnode(np); 
   NgArray<ElementIndex> hasbothpoints;
 
-  BitArray origpoint(np+1), boundp(np+1);  // big enough for 0 and 1-based
+  NgBitArray origpoint(np+1), boundp(np+1);  // big enough for 0 and 1-based
   origpoint.Set();
 
   NgArray<double> elerrs(ne);
-  BitArray illegaltet(ne);
+  NgBitArray illegaltet(ne);
   illegaltet.Clear();
 
   const char * savetask = multithread.task;
@@ -891,7 +891,7 @@ void MeshOptimize3d :: SplitImprove (Mesh & mesh,
   
 
 void MeshOptimize3d :: SwapImprove (Mesh & mesh, OPTIMIZEGOAL goal,
-				    const BitArray * working_elements)
+				    const NgBitArray * working_elements)
 {
   static Timer t("MeshOptimize3d::SwapImprove"); RegionTimer reg(t);
   static Timer tloop("MeshOptimize3d::SwapImprove loop");
@@ -1769,7 +1769,7 @@ void MeshOptimize3d :: SwapImprove (Mesh & mesh, OPTIMIZEGOAL goal,
 
 
 void MeshOptimize3d :: SwapImproveSurface (Mesh & mesh, OPTIMIZEGOAL goal,
-					   const BitArray * working_elements,
+					   const NgBitArray * working_elements,
 					   const NgArray< NgArray<int,PointIndex::BASE>* > * idmaps)
 {
   NgArray< NgArray<int,PointIndex::BASE>* > locidmaps;
@@ -2970,7 +2970,7 @@ void MeshOptimize3d :: SwapImprove2 (Mesh & mesh, OPTIMIZEGOAL goal)
   }
   }
 
-  BitArray original(GetNE());
+  NgBitArray original(GetNE());
   original.Set();
 
   for (i = 1; i <= GetNSE(); i++)

libsrc/meshing/improve3.hpp

@@ -24,9 +24,9 @@ public:
 
   void SplitImprove (Mesh & mesh, OPTIMIZEGOAL goal = OPT_QUALITY);
   void SwapImprove (Mesh & mesh, OPTIMIZEGOAL goal = OPT_QUALITY,
-		    const BitArray * working_elements = NULL);
+		    const NgBitArray * working_elements = NULL);
   void SwapImproveSurface (Mesh & mesh, OPTIMIZEGOAL goal = OPT_QUALITY,
-			   const BitArray * working_elements = NULL,
+			   const NgBitArray * working_elements = NULL,
 			   const NgArray< NgArray<int,PointIndex::BASE>* > * idmaps = NULL);
   void SwapImprove2 (Mesh & mesh, OPTIMIZEGOAL goal = OPT_QUALITY);
 

libsrc/meshing/meshclass.cpp

@@ -1824,8 +1824,8 @@ namespace netgen
       }
   }
 
-  // BitArray base is PointIndex::BASE ... 
+  // NgBitArray base is PointIndex::BASE ... 
-  void Mesh :: FixPoints (const BitArray & fixpoints)
+  void Mesh :: FixPoints (const NgBitArray & fixpoints)
   {
     if (fixpoints.Size() != GetNP())
       {
@@ -2469,7 +2469,7 @@ namespace netgen
     int np = GetNP();
 
     FindOpenSegments();
-    BitArray frontpoints(np+1);  // for 0- and 1-based
+    NgBitArray frontpoints(np+1);  // for 0- and 1-based
     frontpoints.Clear();
     
     for (int i = 1; i <= GetNOpenSegments(); i++)
@@ -2994,7 +2994,7 @@ namespace netgen
     int nse = GetNSE();
 
     NgArray<Vec3d> normals(np);
-    BitArray linepoint(np);
+    NgBitArray linepoint(np);
 
     linepoint.Clear();
     for (i = 1; i <= nseg; i++)
@@ -3260,7 +3260,7 @@ namespace netgen
     
     NgArray<PointIndex,PointIndex::BASE,PointIndex> op2np(GetNP());
     NgArray<MeshPoint> hpoints;
-    BitArrayChar<PointIndex::BASE> pused(GetNP());
+    Array<bool, PointIndex> pused(GetNP());
 
     /*
       (*testout) << "volels: " << endl;
@@ -3300,37 +3300,37 @@ namespace netgen
       if(segments[i].edgenr < 0)
           segments.DeleteElement(i--);
 
-    pused.Clear();
+    pused = false;
     for (int i = 0; i < volelements.Size(); i++)
       {
         const Element & el = volelements[i];
         for (int j = 0; j < el.GetNP(); j++)
-          pused.Set (el[j]);
+          pused[el[j]] = true;
       }
 
     for (int i = 0; i < surfelements.Size(); i++)
       {
         const Element2d & el = surfelements[i];
         for (int j = 0; j < el.GetNP(); j++)
-          pused.Set (el[j]);
+          pused[el[j]] = true;
       }
 
     for (int i = 0; i < segments.Size(); i++)
       {
         const Segment & seg = segments[i];
         for (int j = 0; j < seg.GetNP(); j++)
-          pused.Set (seg[j]);
+          pused[seg[j]] = true;
       }
 
     for (int i = 0; i < openelements.Size(); i++)
       {
         const Element2d & el = openelements[i];
         for (int j = 0; j < el.GetNP(); j++)
-          pused.Set(el[j]);
+          pused[el[j]] = true;
       }
 
     for (int i = 0; i < lockedpoints.Size(); i++)
-      pused.Set (lockedpoints[i]);
+      pused[lockedpoints[i]] = true;
 
 
     /*
@@ -3352,7 +3352,7 @@ namespace netgen
 
     // for (PointIndex pi = points.Begin(); pi < points.End(); pi++)
     for (PointIndex pi : points.Range())
-      if (pused.Test(pi))
+      if (pused[pi])
         {
           npi++;
           op2np[pi] = npi;
@@ -3899,7 +3899,7 @@ namespace netgen
     int i, j;
     int nse = GetNSE();
 
-    BitArray used(nse);
+    NgBitArray used(nse);
     used.Clear();
     INDEX_2_HASHTABLE<int> edges(nse+1);
 
@@ -5289,8 +5289,8 @@ namespace netgen
     int np = GetNP();
     int nse = GetNSE();
 
-    BitArray surfused(nse);
+    NgBitArray surfused(nse);
-    BitArray pused (np);
+    NgBitArray pused (np);
 
     surfused.Clear();
 
@@ -5401,7 +5401,7 @@ namespace netgen
     int fdi;
     int np = GetNP();
 
-    BitArray usedp(np);
+    NgBitArray usedp(np);
     Array<SurfaceElementIndex> els_of_face;
 
     fdi = 1;

libsrc/meshing/meshclass.hpp

@@ -381,7 +381,7 @@ namespace netgen
     DLL_HEADER void CalcSurfacesOfNode ();
 
     /// additional (temporarily) fix points 
-    void FixPoints (const BitArray & fixpoints);
+    void FixPoints (const NgBitArray & fixpoints);
 
     /**
        finds elements without neighbour and
@@ -544,10 +544,10 @@ namespace netgen
 	DLL_HEADER void ImproveMesh (const MeshingParameters & mp, OPTIMIZEGOAL goal = OPT_QUALITY);
 
     ///
-    void ImproveMeshJacobian (const MeshingParameters & mp, OPTIMIZEGOAL goal = OPT_QUALITY, const BitArray * usepoint = NULL);
+    void ImproveMeshJacobian (const MeshingParameters & mp, OPTIMIZEGOAL goal = OPT_QUALITY, const NgBitArray * usepoint = NULL);
     ///
     void ImproveMeshJacobianOnSurface (const MeshingParameters & mp,
-				       const BitArray & usepoint, 
+				       const NgBitArray & usepoint, 
 				       const NgArray< Vec<3>* > & nv,
 				       OPTIMIZEGOAL goal = OPT_QUALITY,
 				       const NgArray< NgArray<int,PointIndex::BASE>* > * idmaps = NULL);

libsrc/meshing/meshtool.cpp

@@ -964,7 +964,7 @@ namespace netgen
     mesh.FindOpenElements(domainnr);
     int np = mesh.GetNP();
 
-    BitArrayChar<PointIndex::BASE> ppoints(np);
+    Array<bool, PointIndex> ppoints(np);
   
     // int ndom = mesh.GetNDomains();
 
@@ -972,7 +972,7 @@ namespace netgen
     // for (k = 1; k <= ndom; k++)
     k = domainnr;
       {
-	ppoints.Clear();
+	ppoints = false;
       
 	for (i = 1; i <= mesh.GetNOpenElements(); i++)
 	  {
@@ -980,7 +980,7 @@ namespace netgen
 	    if (sel.GetIndex() == k)
 	      {
 		for (j = 1; j <= sel.GetNP(); j++)
-		  ppoints.Set (sel.PNum(j));
+		  ppoints[sel.PNum(j)] = true;
 	      }
 	  }
 
@@ -991,7 +991,7 @@ namespace netgen
 	      {
 		int todel = 0;
 		for (j = 0; j < el.GetNP(); j++)
-		  if (ppoints.Test (el[j]))
+		  if (ppoints[el[j]])
 		    todel = 1;
 	      
 		if (el.GetNP() != 4)

libsrc/meshing/refine.cpp

@@ -770,7 +770,7 @@ namespace netgen
 					can.Elem(parent.I2()));
 	    }
 
-	BitArray boundp(np);
+	NgBitArray boundp(np);
 	boundp.Clear();
 	for (auto & sel : mesh.SurfaceElements())
           for (auto pi : sel.PNums())
@@ -801,7 +801,7 @@ namespace netgen
 		    mesh.Point(i) = can.Get(i);
 	      
 
-		BitArray free (mesh.GetNP()), fhelp(mesh.GetNP());
+		NgBitArray free (mesh.GetNP()), fhelp(mesh.GetNP());
 		free.Clear();
 		for (int i = 1; i <= mesh.GetNE(); i++)
 		  {

libsrc/meshing/secondorder.cpp

@@ -480,7 +480,7 @@ namespace netgen
     double facok = 0;
     double factry;
 
-    BitArray illegalels(ne);
+    NgBitArray illegalels(ne);
     illegalels.Clear();
 
       
@@ -504,7 +504,7 @@ namespace netgen
 				    can.Elem(parents.Get(i).I2()));
 	  }
 
-	BitArray boundp(np);
+	NgBitArray boundp(np);
 	boundp.Clear();
 	for (int i = 1; i <= mesh.GetNSE(); i++)
 	  {

libsrc/meshing/smoothing2.5.cpp

@@ -110,7 +110,7 @@ namespace netgen
     int np = mesh.GetNP();
     int ne = mesh.GetNE();
 
-    BitArray badnodes(np);
+    NgBitArray badnodes(np);
     badnodes.Clear();
 
     for (i = 1; i <= ne; i++)

libsrc/meshing/smoothing3.cpp

@@ -1485,7 +1485,7 @@ void Mesh :: ImproveMesh (const MeshingParameters & mp, OPTIMIZEGOAL goal)
 
 // Improve Condition number of Jacobian, any elements  
 void Mesh :: ImproveMeshJacobian (const MeshingParameters & mp,
-				  OPTIMIZEGOAL goal, const BitArray * usepoint)
+				  OPTIMIZEGOAL goal, const NgBitArray * usepoint)
 {
   // int i, j;
   
@@ -1507,7 +1507,7 @@ void Mesh :: ImproveMeshJacobian (const MeshingParameters & mp,
   par.maxit_linsearch = 20;
   par.maxit_bfgs = 20;
   
-  BitArray badnodes(np);
+  NgBitArray badnodes(np);
   badnodes.Clear();
 
   for (int i = 1; i <= ne; i++)
@@ -1608,7 +1608,7 @@ void Mesh :: ImproveMeshJacobian (const MeshingParameters & mp,
 
 // Improve Condition number of Jacobian, any elements  
 void Mesh :: ImproveMeshJacobianOnSurface (const MeshingParameters & mp,
-					   const BitArray & usepoint, 
+					   const NgBitArray & usepoint, 
 					   const NgArray< Vec<3>* > & nv,
 					   OPTIMIZEGOAL goal,
 					   const NgArray< NgArray<int,PointIndex::BASE>* > * idmaps)
@@ -1664,7 +1664,7 @@ void Mesh :: ImproveMeshJacobianOnSurface (const MeshingParameters & mp,
   par.maxit_linsearch = 20;
   par.maxit_bfgs = 20;
   
-  BitArray badnodes(np);
+  NgBitArray badnodes(np);
   badnodes.Clear();
 
   for (int i = 1; i <= ne; i++)

libsrc/meshing/specials.cpp

@@ -62,7 +62,7 @@ void CutOffAndCombine (Mesh & mesh, const Mesh & othermesh)
     }
   cout << endl;
 
-  BitArray connected(mesh.GetNP());
+  NgBitArray connected(mesh.GetNP());
   connected.Clear();
   for (i = 1; i <= mesh.GetNSE(); i++)
     {
@@ -120,7 +120,7 @@ void CutOffAndCombine (Mesh & mesh, const Mesh & othermesh)
   mesh.Compress();
   
   mesh.FindOpenElements();
-  BitArray locked(mesh.GetNP());
+  NgBitArray locked(mesh.GetNP());
   locked.Set();
   for (i = 1; i <= mesh.GetNOpenElements(); i++)
     for (j = 1; j <= 3; j++)

libsrc/meshing/validate.cpp

@@ -6,7 +6,7 @@
 namespace netgen
 {
   void GetPureBadness(Mesh & mesh, NgArray<double> & pure_badness,
-		      const BitArray & isnewpoint)
+		      const NgBitArray & isnewpoint)
   {
     //const int ne = mesh.GetNE();
     const int np = mesh.GetNP();
@@ -104,7 +104,7 @@ namespace netgen
   }
 
 
-  void GetWorkingArea(BitArray & working_elements, BitArray & working_points,
+  void GetWorkingArea(NgBitArray & working_elements, NgBitArray & working_points,
 		      const Mesh & mesh, const NgArray<ElementIndex> & bad_elements,
 		      const int width)
   {
@@ -152,7 +152,7 @@ namespace netgen
 
 
   void RepairBisection(Mesh & mesh, NgArray<ElementIndex> & bad_elements, 
-		       const BitArray & isnewpoint, const Refinement & refinement,
+		       const NgBitArray & isnewpoint, const Refinement & refinement,
 		       const NgArray<double> & pure_badness, 
 		       double max_worsening, const bool uselocalworsening,
 		       const NgArray< NgArray<int,PointIndex::BASE>* > & idmaps)
@@ -185,7 +185,7 @@ namespace netgen
 	can[i] = new Point<3>;
       }
     
-    BitArray isboundarypoint(np),isedgepoint(np);
+    NgBitArray isboundarypoint(np),isedgepoint(np);
     isboundarypoint.Clear();
     isedgepoint.Clear();
 
@@ -216,8 +216,8 @@ namespace netgen
     Validate(mesh,bad_elements,pure_badness,
 	     ((uselocalworsening) ?  (0.8*(max_worsening-1.) + 1.) : (0.1*(max_worsening-1.) + 1.)),
 	     uselocalworsening); // -> larger working area
-    BitArray working_elements(ne);
+    NgBitArray working_elements(ne);
-    BitArray working_points(np);
+    NgBitArray working_points(np);
 
     GetWorkingArea(working_elements,working_points,mesh,bad_elements,numbadneighbours);
     //working_elements.Set();
@@ -240,7 +240,7 @@ namespace netgen
     PrintMessage(5,ostrstr.str());
     
 
-    BitArray isworkingboundary(np);
+    NgBitArray isworkingboundary(np);
     for(int i=1; i<=np; i++)
       if(working_points.Test(i) && isboundarypoint.Test(i))
 	isworkingboundary.Set(i);

libsrc/meshing/validate.hpp

@@ -5,13 +5,13 @@ namespace netgen
 {
   
   void GetPureBadness(Mesh & mesh, NgArray<double> & pure_badness,
-		      const BitArray & isnewpoint);
+		      const NgBitArray & isnewpoint);
   double Validate(const Mesh & mesh, NgArray<ElementIndex> & bad_elements,
 		  const NgArray<double> & pure_badness, 
 		  double max_worsening, const bool uselocalworsening,
 		  NgArray<double> * quality_loss = NULL);
   void RepairBisection(Mesh & mesh, NgArray<ElementIndex> & bad_elements, 
-		       const BitArray & isnewpoint, const Refinement & refinement,
+		       const NgBitArray & isnewpoint, const Refinement & refinement,
 		       const NgArray<double> & pure_badness, 
 		       double max_worsening, const bool uselocalworsening,
 		       const NgArray< NgArray<int,PointIndex::BASE>* > & idmaps);

libsrc/visualization/vsmesh.cpp

@@ -1819,7 +1819,7 @@ namespace netgen
 
     NgArray<Element2d> faces;
 
-    BitArray shownode(mesh->GetNP());
+    NgBitArray shownode(mesh->GetNP());
     if (vispar.clipping.enable)
       {
 	shownode.Clear();

tests/pytest/test_bitarray.py

@@ -0,0 +1,35 @@
+from pyngcore import BitArray
+
+def test_bitarray():
+    a = BitArray(498)
+    assert len(a) == 498
+
+    a.Set()
+    for b in a:
+        assert b == True
+
+    a.Clear(23)
+    assert a[22] == True
+    assert a[23] == False
+    assert a[24] == True
+
+    a.Clear()
+    for b in a:
+        assert b == False
+
+    a.Set(23)
+    assert a[22] == False
+    assert a[23] == True
+    assert a[24] == False
+
+    a.Clear()
+    a[100:200:9] = True
+    for i in range(len(a)):
+        assert a[i] == bool(100<=i and i<200 and i%9==100%9)
+
+    ac = ~a
+
+    for b,bc in zip(a,ac):
+        assert b == (not bc)
+
+