()
| 123 | return pymesh.form_mesh(vertices, tris, tets); |
| 124 | |
| 125 | def main(): |
| 126 | args = parse_args(); |
| 127 | mesh = pymesh.load_mesh(args.input_mesh); |
| 128 | num_tets = mesh.num_voxels; |
| 129 | assert(mesh.dim == 3); |
| 130 | assert(mesh.vertex_per_voxel == 4); |
| 131 | assert(num_tets > 0); |
| 132 | mesh = fit_into_unit_sphere(mesh); |
| 133 | mesh.add_attribute("voxel_centroid"); |
| 134 | assembler = pymesh.Assembler(mesh); |
| 135 | |
| 136 | if args.charge_distribution == "sphere": |
| 137 | ball = pymesh.generate_icosphere(1.25 + 1e-3, np.zeros(3), 0); |
| 138 | charges = ball.vertices; |
| 139 | elif args.charge_distribution == "center": |
| 140 | charges = np.zeros(3); |
| 141 | else: |
| 142 | raise NotImplementedError("Unsupported charge distribution ({})."\ |
| 143 | .format(args.charge_distribution)); |
| 144 | |
| 145 | if args.output_kernels: |
| 146 | kernels = pymesh.merge_meshes([ |
| 147 | pymesh.generate_icosphere(0.05, v, 2) for v in charges ]); |
| 148 | pymesh.save_mesh("kernel.msh", kernels); |
| 149 | f = test_function(charges); |
| 150 | g = test_function_grad(charges); |
| 151 | |
| 152 | v_pts = mesh.vertices; |
| 153 | q_pts = get_quadrature_pts(mesh); |
| 154 | c_pts = mesh.get_voxel_attribute("voxel_centroid"); |
| 155 | |
| 156 | v_values = f(v_pts); |
| 157 | q_values = f(q_pts); |
| 158 | c_values = f(c_pts); |
| 159 | q_grad = g(q_pts); |
| 160 | c_grad = g(c_pts); |
| 161 | c_grad_norm = numpy.linalg.norm(c_grad, axis=1); |
| 162 | q_grad_norm = numpy.linalg.norm(q_grad, axis=1); |
| 163 | #for v,val in zip(v_pts, v_values): |
| 164 | # print(v, val); |
| 165 | |
| 166 | assert(len(v_values) == len(v_pts)); |
| 167 | assert(len(q_values) == len(q_pts)); |
| 168 | assert(len(c_values) == len(c_pts)); |
| 169 | assert(len(q_grad) == len(q_pts)); |
| 170 | assert(len(c_grad) == len(c_pts)); |
| 171 | |
| 172 | bd_indices = np.unique(mesh.faces.ravel()).astype(int); |
| 173 | bd_values = v_values[bd_indices]; |
| 174 | |
| 175 | solver = pymesh.HarmonicSolver.create(mesh); |
| 176 | solver.order = 1; |
| 177 | solver.boundary_indices = bd_indices; |
| 178 | solver.boundary_values = bd_values; |
| 179 | solver.pre_process(); |
| 180 | solver.solve(); |
| 181 | |
| 182 | sol_values = solver.solution.ravel(); |
no test coverage detected