1
1
import numpy as np
2
- from discretize import utils
2
+ import properties
3
+ import os
4
+ import json
5
+ from .utils .matutils import mkvc
3
6
4
7
5
- class BaseMesh (object ):
6
- """BaseMesh does all the counting you don't want to do.
8
+ class BaseMesh (properties .HasProperties ):
9
+ """
10
+ BaseMesh does all the counting you don't want to do.
7
11
BaseMesh should be inherited by meshes with a regular structure.
8
-
9
- :param numpy.array n: (or list) number of cells in each direction (dim, )
10
- :param numpy.array x0: (or list) Origin of the mesh (dim, )
11
-
12
12
"""
13
13
14
- def __init__ ( self , n , x0 = None ):
14
+ _REGISTRY = {}
15
15
16
- # Check inputs
17
- if x0 is None :
18
- x0 = np .zeros (len (n ))
16
+ # Properties
17
+ _n = properties .Array (
18
+ "number of cells in each direction (dim, )" ,
19
+ dtype = int ,
20
+ required = True ,
21
+ shape = ('*' ,)
22
+ )
19
23
20
- if not len (n ) == len (x0 ):
21
- raise Exception ("Dimension mismatch. x0 != len(n)" )
24
+ x0 = properties .Array (
25
+ "origin of the mesh (dim, )" ,
26
+ dtype = float ,
27
+ shape = ('*' ,),
28
+ required = True
29
+ )
22
30
23
- if len (n ) > 3 :
24
- raise Exception ("Dimensions higher than 3 are not supported." )
31
+ # Instantiate the class
32
+ def __init__ (self , n , x0 = None ):
33
+ self ._n = n # number of dimensions
25
34
26
- # Ensure x0 & n are 1D vectors
27
- self ._n = np .array ( n , dtype = int ). ravel ( )
28
- self . _x0 = np . array ( x0 , dtype = float ). ravel ()
29
- self ._dim = len ( self . _x0 )
35
+ if x0 is None :
36
+ self .x0 = np .zeros ( len ( self . _n ) )
37
+ else :
38
+ self .x0 = x0
30
39
31
- @property
32
- def x0 (self ):
33
- """Origin of the mesh
40
+ # Validators
41
+ @properties .validator ('_n' )
42
+ def check_n_shape (self , change ):
43
+ assert (
44
+ not isinstance (change ['value' ], properties .utils .Sentinel ) and
45
+ change ['value' ] is not None
46
+ ), "Cannot delete n. Instead, create a new mesh"
47
+
48
+ change ['value' ] = np .array (change ['value' ], dtype = int ).ravel ()
49
+ if len (change ['value' ]) > 3 :
50
+ raise Exception (
51
+ "Dimensions of {}, which is higher than 3 are not "
52
+ "supported" .format (change ['value' ])
53
+ )
54
+
55
+ if change ['previous' ] != properties .undefined :
56
+ # can't change dimension of the mesh
57
+ assert len (change ['previous' ]) == len (change ['value' ]), (
58
+ "Cannot change dimensionality of the mesh. Expected {} "
59
+ "dimensions, got {} dimensions" .format (
60
+ len (change ['previous' ]), len (change ['value' ])
61
+ )
62
+ )
63
+
64
+ # check that if h has been set, sizes still agree
65
+ if getattr (self , 'h' , None ) is not None and len (self .h ) > 0 :
66
+ for i in range (len (change ['value' ])):
67
+ assert len (self .h [i ]) == change ['value' ][i ], (
68
+ "Mismatched shape of n. Expected {}, len(h[{}]), got "
69
+ "{}" .format (
70
+ len (self .h [i ]), i , change ['value' ][i ]
71
+ )
72
+ )
73
+
74
+ # check that if nodes have been set for curvi mesh, sizes still
75
+ # agree
76
+ if (
77
+ getattr (self , 'nodes' , None ) is not None and
78
+ len (self .nodes ) > 0
79
+ ):
80
+ for i in range (len (change ['value' ])):
81
+ assert self .nodes [0 ].shape [i ]- 1 == change ['value' ][i ], (
82
+ "Mismatched shape of n. Expected {}, len(nodes[{}]), "
83
+ "got {}" .format (
84
+ self .nodes [0 ].shape [i ]- 1 , i , change ['value' ][i ]
85
+ )
86
+ )
87
+
88
+ @properties .validator ('x0' )
89
+ def check_x0 (self , change ):
90
+ assert (
91
+ not isinstance (change ['value' ], properties .utils .Sentinel ) and
92
+ change ['value' ] is not None
93
+ ), "n must be set prior to setting x0"
34
94
35
- :rtype: numpy.array
36
- :return: x0, (dim, )
37
- """
38
- return self ._x0
95
+ if len (self ._n ) != len (change ['value' ]):
96
+ raise Exception (
97
+ "Dimension mismatch. x0 has length {} != len(n) which is "
98
+ "{}" .format (len (x0 ), len (n ))
99
+ )
39
100
40
101
@property
41
102
def dim (self ):
@@ -44,7 +105,7 @@ def dim(self):
44
105
:rtype: int
45
106
:return: dim
46
107
"""
47
- return self ._dim
108
+ return len ( self ._n )
48
109
49
110
@property
50
111
def nC (self ):
@@ -290,11 +351,33 @@ def projectEdgeVector(self, eV):
290
351
), 'eV must be an ndarray of shape (nE x dim)'
291
352
return np .sum (eV * self .tangents , 1 )
292
353
354
+ def save (self , filename = 'mesh.json' , verbose = False ):
355
+ """
356
+ Save the mesh to json
357
+ :param str file: filename for saving the casing properties
358
+ :param str directory: working directory for saving the file
359
+ """
360
+
361
+ f = os .path .abspath (filename ) # make sure we are working with abs path
362
+ with open (f , 'w' ) as outfile :
363
+ json .dump (self .serialize (), outfile )
364
+
365
+ if verbose is True :
366
+ print ('Saved {}' .format (f ))
367
+
368
+ return f
369
+
370
+ def copy (self ):
371
+ """
372
+ Make a copy of the current mesh
373
+ """
374
+ return properties .copy (self )
375
+
293
376
294
377
class BaseRectangularMesh (BaseMesh ):
295
378
"""BaseRectangularMesh"""
296
379
def __init__ (self , n , x0 = None ):
297
- BaseMesh .__init__ (self , n , x0 )
380
+ BaseMesh .__init__ (self , n , x0 = x0 )
298
381
299
382
@property
300
383
def nCx (self ):
@@ -590,44 +673,69 @@ def r(self, x, xType='CC', outType='CC', format='V'):
590
673
eX, eY, eZ = r(edgeVector, 'E', 'E', 'V')
591
674
"""
592
675
593
- assert (type (x ) == list or isinstance (x , np .ndarray )), "x must be either a list or a ndarray"
594
- assert xType in ['CC' , 'N' , 'F' , 'Fx' , 'Fy' , 'Fz' , 'E' , 'Ex' , 'Ey' , 'Ez' ], "xType must be either 'CC', 'N', 'F', 'Fx', 'Fy', 'Fz', 'E', 'Ex', 'Ey', or 'Ez'"
595
- assert outType in ['CC' , 'N' , 'F' , 'Fx' , 'Fy' , 'Fz' , 'E' , 'Ex' , 'Ey' , 'Ez' ], "outType must be either 'CC', 'N', 'F', Fx', 'Fy', 'Fz', 'E', 'Ex', 'Ey', or 'Ez'"
676
+ allowed_xType = [
677
+ 'CC' , 'N' , 'F' , 'Fx' , 'Fy' , 'Fz' , 'E' , 'Ex' , 'Ey' , 'Ez'
678
+ ]
679
+ assert (
680
+ type (x ) == list or isinstance (x , np .ndarray )
681
+ ), "x must be either a list or a ndarray"
682
+ assert xType in allowed_xType , (
683
+ "xType must be either "
684
+ "'CC', 'N', 'F', 'Fx', 'Fy', 'Fz', 'E', 'Ex', 'Ey', or 'Ez'"
685
+ )
686
+ assert outType in allowed_xType , (
687
+ "outType must be either "
688
+ "'CC', 'N', 'F', Fx', 'Fy', 'Fz', 'E', 'Ex', 'Ey', or 'Ez'"
689
+ )
596
690
assert format in ['M' , 'V' ], "format must be either 'M' or 'V'"
597
- assert outType [:len (xType )] == xType , "You cannot change types when reshaping."
598
- assert xType in outType , 'You cannot change type of components.'
691
+ assert outType [:len (xType )] == xType , (
692
+ "You cannot change types when reshaping."
693
+ )
694
+ assert xType in outType , "You cannot change type of components."
695
+
599
696
if type (x ) == list :
600
697
for i , xi in enumerate (x ):
601
- assert isinstance (x , np .ndarray ), "x[{0:d}] must be a numpy array" .format (i )
602
- assert xi .size == x [0 ].size , "Number of elements in list must not change."
698
+ assert isinstance (x , np .ndarray ), (
699
+ "x[{0:d}] must be a numpy array" .format (i )
700
+ )
701
+ assert xi .size == x [0 ].size , (
702
+ "Number of elements in list must not change."
703
+ )
603
704
604
705
x_array = np .ones ((x .size , len (x )))
605
706
# Unwrap it and put it in a np array
606
707
for i , xi in enumerate (x ):
607
- x_array [:, i ] = utils . mkvc (xi )
708
+ x_array [:, i ] = mkvc (xi )
608
709
x = x_array
609
710
610
711
assert isinstance (x , np .ndarray ), "x must be a numpy array"
611
712
612
713
x = x [:] # make a copy.
613
- xTypeIsFExyz = len (xType ) > 1 and xType [0 ] in ['F' , 'E' ] and xType [1 ] in ['x' , 'y' , 'z' ]
714
+ xTypeIsFExyz = (
715
+ len (xType ) > 1 and
716
+ xType [0 ] in ['F' , 'E' ] and
717
+ xType [1 ] in ['x' , 'y' , 'z' ]
718
+ )
614
719
615
720
def outKernal (xx , nn ):
616
721
"""Returns xx as either a matrix (shape == nn) or a vector."""
617
722
if format == 'M' :
618
723
return xx .reshape (nn , order = 'F' )
619
724
elif format == 'V' :
620
- return utils . mkvc (xx )
725
+ return mkvc (xx )
621
726
622
727
def switchKernal (xx ):
623
728
"""Switches over the different options."""
624
729
if xType in ['CC' , 'N' ]:
625
730
nn = (self ._n ) if xType == 'CC' else (self ._n + 1 )
626
- assert xx .size == np .prod (nn ), "Number of elements must not change."
731
+ assert xx .size == np .prod (nn ), (
732
+ "Number of elements must not change."
733
+ )
627
734
return outKernal (xx , nn )
628
735
elif xType in ['F' , 'E' ]:
629
- # This will only deal with components of fields, not full 'F' or 'E'
630
- xx = utils .mkvc (xx ) # unwrap it in case it is a matrix
736
+ # This will only deal with components of fields,
737
+ # not full 'F' or 'E'
738
+ xx = mkvc (xx ) # unwrap it in case it is a matrix
631
739
nn = self .vnF if xType == 'F' else self .vnE
632
740
nn = np .r_ [0 , nn ]
633
741
@@ -638,13 +746,20 @@ def switchKernal(xx):
638
746
639
747
for dim , dimName in enumerate (['x' , 'y' , 'z' ]):
640
748
if dimName in outType :
641
- assert self .dim > dim , ("Dimensions of mesh not great enough for %s%s" , (xType , dimName ))
642
- assert xx .size == np .sum (nn ), 'Vector is not the right size.'
749
+ assert self .dim > dim , (
750
+ "Dimensions of mesh not great enough for "
751
+ "{}{}" .format (xType , dimName )
752
+ )
753
+ assert xx .size == np .sum (nn ), (
754
+ "Vector is not the right size."
755
+ )
643
756
start = np .sum (nn [:dim + 1 ])
644
757
end = np .sum (nn [:dim + 2 ])
645
758
return outKernal (xx [start :end ], nx [dim ])
759
+
646
760
elif xTypeIsFExyz :
647
- # This will deal with partial components (x, y or z) lying on edges or faces
761
+ # This will deal with partial components (x, y or z)
762
+ # lying on edges or faces
648
763
if 'x' in xType :
649
764
nn = self .vnFx if 'F' in xType else self .vnEx
650
765
elif 'y' in xType :
@@ -658,7 +773,9 @@ def switchKernal(xx):
658
773
isVectorQuantity = len (x .shape ) == 2 and x .shape [1 ] == self .dim
659
774
660
775
if outType in ['F' , 'E' ]:
661
- assert ~ isVectorQuantity , 'Not sure what to do with a vector vector quantity..'
776
+ assert ~ isVectorQuantity , (
777
+ 'Not sure what to do with a vector vector quantity..'
778
+ )
662
779
outTypeCopy = outType
663
780
out = ()
664
781
for ii , dirName in enumerate (['x' , 'y' , 'z' ][:self .dim ]):
0 commit comments