Matthias/refactoring rkpm weak dirichlet locations

.github/workflows/run_tests.yml

@@ -10,13 +10,13 @@ jobs:
       - uses: actions/checkout@v2
 
       - name: Setup Mambaforge
-        uses: conda-incubator/setup-miniconda@v2
+        uses: conda-incubator/setup-miniconda@v3
         with:
           use-mamba: true
-          miniforge-variant: Mambaforge
+          miniforge-variant: Miniforge3
           miniforge-version: latest
           activate-environment: anaconda-client-env
-              
+
       - name: Install conda packages
         shell: bash -l {0}
         run: |
@@ -27,9 +27,9 @@ jobs:
         shell: bash -l {0}
         run: |
           git clone --branch 3.4.0 https://gitlab.com/libeigen/eigen.git
-          cd eigen 
+          cd eigen
           mkdir build
-          cd build 
+          cd build
           cmake -DBUILD_TESTING=OFF -DCMAKE_INSTALL_PREFIX=$CONDA_PREFIX -DINCLUDE_INSTALL_DIR=$CONDA_PREFIX/include ..
           make install
           cd ../..
@@ -57,16 +57,13 @@ jobs:
       - name: Install Marmot
         shell: bash -l {0}
         run: |
-          git clone --recurse-submodules https://github.com/MAteRialMOdelingToolbox/Marmot/
-          cd Marmot
-          cd modules/core
-          git clone ${{secrets.MARMOT_MPMCORE_TOKEN}}
-          cd ../../
-          mkdir build
-          cd build
+          git clone "https://github.com/MAteRialMOdelingToolbox/Marmot.git"
+          git clone "https://github.com/BOKU-CMP-MEC-MAT/MarmotMeshfreeCore.git" ./Marmot/modules/core/MarmotMeshfreeCore
+          mkdir ./Marmot/build
+          cd ./Marmot/build
           cmake -DCMAKE_INSTALL_PREFIX=$CONDA_PREFIX  ..
           make install
-          cd ../../
+          cd ../..
 
       - name: Install EdelweissFE
         shell: bash -l {0}

README.md

@@ -1,5 +1,5 @@
 <!-- [![documentation](https://github.com/EdelweissMPM/EdelweissMPM/actions/workflows/sphinx.yml/badge.svg)](https://edelweissfe.github.io/EdelweissMPM) -->
-[![Code style: black](https://img.shields.io/badge/code%20style-black-000000.svg)](https://github.com/psf/black) 
+[![Code style: black](https://img.shields.io/badge/code%20style-black-000000.svg)](https://github.com/psf/black)
 
 # EdelweissMPM: A light-weight, platform-independent, parallel material point module for EdelweissFE.
 
@@ -14,7 +14,7 @@ Some features are:
 
  * Python for non performance-critical routines
  * Cython for performance-critical routines
- * Parallelization 
+ * Parallelization
  * Modular system, which is easy to extend
  * Output to Paraview, Ensight, CSV, matplotlib
  * Interfaces to powerful direct and iterative linear solvers

doc/source/documentation/material_point_managers.rst → doc/source/documentation/mpm/material_point_managers.rst

@@ -18,5 +18,5 @@ Smart Material Point Manager
 .. automodule:: edelweissmpm.mpmmanagers.smartmpmmanager
    :private-members:
    :members:
-   
-   
+
+

doc/source/index.rst

@@ -9,7 +9,7 @@ EdelweissMPM aims to be...
 
  - ... a development and learning environment for constitutive models and the material point method,
  - ... an easy to use tool for coupled problems,
- - ... a very flexible tool for implementing and employing special techniques (e.g., the indirect displacement control technique), 
+ - ... a very flexible tool for implementing and employing special techniques (e.g., the indirect displacement control technique),
    which are often more difficult and time consuming to implement in mature, MPI-parallelized codes,
  - ... an efficient tool for nonlinear simulations up to medium sized problems (~ :math:`10^5` degrees of freedom).
 

edelweissmpm/cells/base/cell.py

@@ -30,7 +30,7 @@
 to the material points. Cells are used to compute the residual and stiffness
 matrices for the global system of equations.
 
-Implementing your own cells can be done easily by subclassing from 
+Implementing your own cells can be done easily by subclassing from
 the abstract base class :class:`~CellBase`.
 """
 

edelweissmpm/cells/marmotcell/bsplinemarmotcell.pyx

@@ -90,7 +90,7 @@ cdef class BSplineMarmotCellWrapper(MarmotCellWrapper):
                                                                     &knotVectorsView[0,0],
                                                                     knotVectorsView.size)
 
-        except IndexError:
+        except ValueError:
             raise NotImplementedError("Marmot cell {:} not found in library.".format(cellType))
 
     def __dealloc__(self):

edelweissmpm/cells/marmotcell/lagrangianmarmotcell.pyx

@@ -87,7 +87,7 @@ cdef class LagrangianMarmotCellWrapper(MarmotCellWrapper):
 
         try:
             self._marmotCell = MarmotCellFactory.createCell( cellType.encode('utf-8'), self._cellNumber, &self._nodeCoordinates[0,0], self._nodeCoordinates.size)
-        except IndexError:
+        except ValueError:
             raise NotImplementedError("Marmot cell {:} not found in library.".format(cellType))
 
     def __dealloc__(self):

edelweissmpm/cells/test/cell.py

@@ -25,7 +25,7 @@
 #  the top level directory of EdelweissMPM.
 #  ---------------------------------------------------------------------
 """
-Implementing your own cells can be done easily by subclassing from 
+Implementing your own cells can be done easily by subclassing from
 the abstract base class :class:`~CellBase`.
 """
 

edelweissmpm/constraints/particlelagrangianweakdirichlet.py

@@ -34,24 +34,35 @@ def __init__(
         self,
         name: str,
         constrainedParticle: BaseParticle,
-        constrainedLocation: int | str,
         field: str,
         prescribedStepDelta: dict,
         model,
+        location: str = "center",
+        faceID: int = None,
+        vertexID: int = None,
     ):
         self._name = name
         self._constrainedParticle = constrainedParticle
         self._field = field
         self._prescribedStepDelta = prescribedStepDelta
         self._fieldSize = getFieldSize(self._field, model.domainSize)
         self._nodes = dict()
-        self._constrainedLocation = constrainedLocation
 
-        if isinstance(constrainedLocation, str):
-            if constrainedLocation != "center":
-                raise ValueError("Constrain must be 'center' or a vertex index.")
-        # elif isinstance(constrainedLocation, int):
-        #     raise ValueError(f"Constrain must be 'center' or a vertex index, got {type(constrainedLocation)}.")
+        def _getConstraintLocation():
+            particle = self._constrainedParticle
+            if location == "center":
+                constraintLocation = particle.getCenterCoordinates()
+            elif location == "face":
+                if faceID is None:
+                    raise ValueError("faceID must be specified when location is 'face'.")
+                constraintLocation = particle.getFaceCoordinates(faceID)
+            elif location == "vertex":
+                if vertexID is None:
+                    raise ValueError("vertexID must be specified when location is 'vertex'.")
+                constraintLocation = particle.getVertexCoordinates(vertexID)
+            return constraintLocation
+
+        self._getConstraintLocation = _getConstraintLocation
 
         self._nLagrangianMultipliers = len(self._prescribedStepDelta)
         self.reactionForce = np.zeros(self._fieldSize)
@@ -119,10 +130,12 @@ def applyConstraint(self, dU_: np.ndarray, PExt: np.ndarray, V: np.ndarray, time
         self.reactionForce.fill(0.0)
         p = self._constrainedParticle
 
-        if self._constrainedLocation == "center":
-            constrainedCoordinates = p.getCenterCoordinates()
-        elif isinstance(self._constrainedLocation, int):
-            constrainedCoordinates = p.getVertexCoordinates()[self._constrainedLocation]
+        # if self._constrainedLocation == "center":
+        #     constrainedCoordinates = p.getCenterCoordinates()
+        # elif isinstance(self._constrainedLocation, int):
+        #     constrainedCoordinates = p.getVertexCoordinates()[self._constrainedLocation]
+
+        constrainedCoordinates = self._getConstraintLocation()
 
         for lag, (i, prescribedComponent) in enumerate(self._prescribedStepDelta.items()):
             # lag is the lagrange multiplier index
@@ -158,15 +171,19 @@ def applyConstraint(self, dU_: np.ndarray, PExt: np.ndarray, V: np.ndarray, time
 def ParticleLagrangianWeakDirichletOnParticleSetFactory(
     baseName: str,
     particleSet: list[BaseParticle],
-    constrainedLocation: int | str,
     field: str,
     prescribedStepDelta: dict,
     model: MPMModel,
+    location: str = "center",
+    faceID: int = None,
+    vertexID: int = None,
 ):
     constraints = dict()
     for i, p in enumerate(particleSet):
         name = f"{baseName}_{i}"
-        constraint = ParticleLagrangianWeakDirichlet(name, p, constrainedLocation, field, prescribedStepDelta, model)
+        constraint = ParticleLagrangianWeakDirichlet(
+            name, p, field, prescribedStepDelta, model, location, faceID, vertexID
+        )
         constraints[name] = constraint
 
     return constraints

edelweissmpm/constraints/particlepenaltyweakdirichtlet.py

@@ -41,7 +41,7 @@ def __init__(
         prescribedStepDelta: dict,
         penaltyParameter: float,
         constrain: str | list[int] = "center",
-        **kwargs
+        **kwargs,
     ):
         self._name = name
         self._model = model

edelweissmpm/fieldoutput/mpresultcollector.pyx

@@ -1,32 +1,30 @@
-#!/usr/bin/env python3
 # -*- coding: utf-8 -*-
 #  ---------------------------------------------------------------------
 #
-#  _____    _      _              _         _____ _____ 
-# | ____|__| | ___| |_      _____(_)___ ___|  ___| ____|
-# |  _| / _` |/ _ \ \ \ /\ / / _ \ / __/ __| |_  |  _|  
-# | |__| (_| |  __/ |\ V  V /  __/ \__ \__ \  _| | |___ 
-# |_____\__,_|\___|_| \_/\_/ \___|_|___/___/_|   |_____|
-#                                                       
-# 
+#  _____    _      _              _         __  __ ____  __  __
+# | ____|__| | ___| |_      _____(_)___ ___|  \/  |  _ \|  \/  |
+# |  _| / _` |/ _ \ \ \ /\ / / _ \ / __/ __| |\/| | |_) | |\/| |
+# | |__| (_| |  __/ |\ V  V /  __/ \__ \__ \ |  | |  __/| |  | |
+# |_____\__,_|\___|_| \_/\_/ \___|_|___/___/_|  |_|_|   |_|  |_|
+#
+#
 #  Unit of Strength of Materials and Structural Analysis
 #  University of Innsbruck,
-#  2017 - today
-# 
+#  2023 - today
+#
 #  Matthias Neuner matthias.neuner@uibk.ac.at
-# 
-#  This file is part of EdelweissFE.
-# 
+#
+#  This file is part of EdelweissMPM.
+#
 #  This library is free software; you can redistribute it and/or
 #  modify it under the terms of the GNU Lesser General Public
 #  License as published by the Free Software Foundation; either
 #  version 2.1 of the License, or (at your option) any later version.
-# 
+#
 #  The full text of the license can be found in the file LICENSE.md at
-#  the top level directory of EdelweissFE.
+#  the top level directory of EdelweissMPM.
 #  ---------------------------------------------------------------------
 
-
 cimport numpy as np
 
 import numpy as np
@@ -35,30 +33,30 @@ from libc.stdlib cimport free, malloc
 
 
 cdef class MaterialPointResultCollector:
-    
+
     cdef public  resultsTable
-    
-    cdef int nMPs, nSize 
+
+    cdef int nMPs, nSize
     cdef double[:, ::1] res_
     cdef double** resultPointers
-    
+
     def __init__(self, materialPoints:list, result:str):
         """
-        A cdef class for collecting materialPoint results (by using the permanent results pointer (i.e., a numpy array) 
+        A cdef class for collecting materialPoint results (by using the permanent results pointer (i.e., a numpy array)
         in large array of all materialPoints and all quadrature points.
 
         Collecting materialPointal results may be a performance critical part.
         This cdef class allows for the efficient gathering.
         A 3D array is assembled if multiple quadrature points are requested (shape ``[materialPoints, quadraturePoints, resultVector]`` )
         or a 2D array for one quadrature point ( shape ``[materialPoints, resultVector]`` ).
-        
+
         Method :func:`~edelweissfe.utils.materialPointresultcollector.MaterialPointResultCollector.getCurrentResults` updates the assembly array and passes it back.
-        
+
         The caller is responsible to make a copy of it, if persistent results are needed!
 
         Parameters
         ----------
-        materialPoints 
+        materialPoints
             The list of materialPoints for which the results should be collected.
         result
             The name of the requested result.
@@ -67,17 +65,17 @@ cdef class MaterialPointResultCollector:
         # hotfix for cython compile error associated with 'range' typing of argument quadraturePoints
         # this is due to a bug in in cython 0.29.xx and should be fixed in cython 3.x.x
         # https://github.com/cython/cython/issues/4002
-        
+
         self.nMPs = len(materialPoints)
 
         # assemble a 2d list of all permanent result arrays (=continously updated np arrays)
         resultsPointerList = [  el.getResultArray(result, getPersistentView=True)  for el in materialPoints ]
         self.nSize = resultsPointerList[0].shape[0]
 
-        
+
         # allocate an equivalent 2D C-array for the pointers to each materialPoints results
         self.resultPointers = <double**> malloc ( sizeof(double*) * self.nMPs )
-        
+
         cdef double* ptr
         cdef double[::1] res
         # fill the 2D C-array of pointers by accessing the materialPoints' resultArrays memoryviews
@@ -86,13 +84,13 @@ cdef class MaterialPointResultCollector:
             res = el
             ptr = <double*> &res[0]
             self.resultPointers[ i  ] = ptr
-        
+
         # initialize the large assembly array
         self.resultsTable = np.empty([self.nMPs, self.nSize]  )
-        
+
         # an internal use only memoryview is created for accesing the assembly
         self.res_ = self.resultsTable
-    
+
     def update(self, ):
         """Update all results."""
 
@@ -102,7 +100,7 @@ cdef class MaterialPointResultCollector:
                     # most inner loop: could also be handled by copying the complete vector at once,
                     # but this version turned out to be faster!
                     self.res_[i,k] = self.resultPointers[i][k]
-    
+
     def getCurrentResults(self,) -> np.ndarray:
         """Update and get current results.
 
@@ -114,7 +112,7 @@ cdef class MaterialPointResultCollector:
 
         self.update()
         return self.resultsTable
-    
+
     def __dealloc__(self):
         free ( self.resultPointers )