Gsolv direct

MartinAchondo · Jul 1, 2024 · 29b55ba · 29b55ba
1 parent f0b645f
commit 29b55ba
Show file tree

Hide file tree

Showing 2 changed files with 11 additions and 3 deletions.
diff --git a/xppbe/Model/Equations.py b/xppbe/Model/Equations.py
@@ -1,3 +1,4 @@
+import numpy as np
 import tensorflow as tf
 
 from .PDE_Model import PBE
@@ -47,8 +48,15 @@ def get_dphi(self,X,Nv,flag,model,value='phi'):
         return du_1,du_2
 
     def get_solvation_energy(self,model):
-        X = tf.reshape(tf.constant(self.grid.vertices.transpose(),dtype=self.DTYPE), (-1,3))
-        Nv = tf.reshape(tf.constant(self.grid.normals.transpose(),dtype=self.DTYPE), (-1,3))
+        vertices = self.grid.vertices
+        faces_normals = self.grid.normals
+        elements = self.grid.elements
+        centroids = np.zeros((3, elements.shape[1]))
+        for i, element in enumerate(elements.T):
+            centroids[:, i] = np.mean(vertices[:, element], axis=1)
+
+        X = tf.reshape(tf.constant(centroids.transpose(),dtype=self.DTYPE), (-1,3))
+        Nv = tf.reshape(tf.constant(faces_normals.transpose(),dtype=self.DTYPE), (-1,3))
         phi = model(X,'interface')
         phi_mean = (phi[:,0]+phi[:,1])/2
 

diff --git a/xppbe/Model/PDE_Model.py b/xppbe/Model/PDE_Model.py
@@ -309,7 +309,7 @@ def get_integral_operators(self):
         elements = self.mesh.mol_faces
         vertices = self.mesh.mol_verts
         self.grid = self.bempp.Grid(vertices.transpose(), elements.transpose())
-        self.space = self.bempp.function_space(self.grid, "P", 1)
+        self.space = self.bempp.function_space(self.grid, "DP", 0)
         self.dirichl_space = self.space
         self.neumann_space = self.space