mumpslsys.py 2.96 KB
Newer Older
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
# MigFlow - Copyright (C) <2010-2018>
# <Universite catholique de Louvain (UCL), Belgium
#  Universite de Montpellier, France>
# 	
# List of the contributors to the development of MigFlow: see AUTHORS file.
# Description and complete License: see LICENSE file.
# 	
# This program (MigFlow) 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
# 3 of the License, or (at your option) any later version.
# 
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
# GNU Lesser General Public License for more details.
# 
# You should have received a copy of the GNU Lesser General Public License
# along with this program (see COPYING and COPYING.LESSER files).  If not, 
# see self.mat<http://www.gnu.org/licenses/>.
from gmsh import option
import numpy as np


import ctypes
import ctypes.util
import os

dir_path = os.path.dirname(os.path.realpath(__file__))
lib2 = np.ctypeslib.load_library("libmbfluid2",dir_path)

options_key = {"verbosity":4, "fill_in":14, "reorder_metis":7, "scaling":8, "openmp_threads" :16}

def _np2c(a,dtype=float,order="C") :
    tmp = np.require(a,dtype,order)
    r = ctypes.c_void_p(tmp.ctypes.data)
    r.tmp = tmp
    return r

class LinearSystem :
    def __init__(self, elements, n_fields, options=None) :
        n_nodes = np.max(elements)+1
        self.globalsize = n_nodes*n_fields
        self.localsize = elements.shape[1]*n_fields
        self.idx = (elements[:,None,:]*n_fields+np.arange(n_fields)[None,:,None]).reshape([-1])
        idxf = ((elements[:,:,None]*n_fields+np.arange(n_fields)[None,None,:]).reshape([-1])+1).astype(np.int32)
        eltptr = np.arange(1,elements.size*n_fields+2,self.localsize)
        lib2.mumps_element_new.restype = ctypes.c_void_p
        self.mumps_ptr = ctypes.c_void_p(lib2.mumps_element_new(ctypes.c_int(self.globalsize),ctypes.c_int(eltptr.size-1),_np2c(eltptr,np.int32),_np2c(idxf,np.int32)))  
        options = None if options=="" else options
        if options is not None:
            options_vec = []
            for key, val in options.items():
                if isinstance(key,str):
                    key = options_key[key]
                    options_vec.extend([key,val])
                    print(options_vec)
            options_vec = np.array(options_vec, dtype=np.int32)
            lib2.mumps_set_options(self.mumps_ptr, _np2c(options_vec, np.int32), ctypes.c_int(np.size(options_vec)//2))


    def local_to_global(self,localv,localm,rhs):
        rhs[:] = np.bincount(self.idx,localv,self.globalsize)
        self.v = localm.reshape([-1])

    def solve(self,rhs) :
        r = np.copy(rhs)
        lib2.mumps_element_solve(self.mumps_ptr,_np2c(self.v),_np2c(r))
        return r

    def __delete__(self):
        lib2.mumps_element_delete(self.mumps_ptr)