1""" 2Create header files with: 3 4- triangle quadrature rules from the supplementary data in https://doi.org/10.1016/j.camwa.2015.03.017 5- tetrahedral quadrature rules from the supplementary data in https://doi.org/10.1002/nme.6528 6""" 7 8import os 9import re 10import glob 11import gmpy2 12from gmpy2 import mpfr 13import pprint 14 15gmpy2.get_context().precision = 113 # Generate quad precision headers 16 17def triangle_rules(rules_path): 18 """ 19 Gather and format quadrature rule data from the unpacked zip file found 20 at https://doi.org/10.1016/j.camwa.2015.03.017 21 """ 22 23 compact_folder = os.path.join(rules_path, "compact", "tri") 24 expanded_folder = os.path.join(rules_path, "expanded", "tri") 25 26 def biunit_to_bary(x): 27 b = [mpfr(0), mpfr(0), mpfr(0)] 28 btob = [ 29 [ 0.5, 0.0, 0.5], 30 [ 0.0, 0.5, 0.5], 31 [-0.5,-0.5, 0.0] 32 ] 33 for i in range(3): 34 for j in range(3): 35 b[i] += btob[i][j] * x[j] 36 return b 37 38 rules = {} 39 for file in os.listdir(compact_folder): 40 name = os.path.basename(file) 41 descrip = os.path.splitext(name)[0] 42 parts = descrip.split('-') 43 power = int(parts[0]) 44 num_points = int(parts[1]) 45 counts = () 46 with open(os.path.join(compact_folder,file),"r") as f: 47 counts = tuple(int(s) for s in re.split('\s+',f.readline().strip())) 48 rule = {} 49 rule['num_nodes'] = num_points 50 rule['counts'] = counts 51 rule['nodes'] = [[], [], []] 52 rule['weights'] = [] 53 with open(os.path.join(expanded_folder,file),"r") as f: 54 for (i,(count, multiplicity)) in enumerate(zip(counts, (1,3,6))): 55 for c in range(count): 56 for m in range(multiplicity): 57 x, y, weight = tuple(mpfr(s) for s in re.split('\s+',f.readline().strip())) 58 if m == 0: 59 rule['weights'].append(weight) 60 if i == 0: 61 rule['nodes'][i].append([mpfr(1) / mpfr(3)]) 62 elif i == 1: 63 X = [x, y, mpfr(1)] 64 B = biunit_to_bary(X) 65 nodes = sorted(B) 66 diffs = [b - a for (a,b) in zip(nodes[:-1],nodes[1:])] 67 if diffs[0] > diffs[1]: # the first is the odd one out 68 rule['nodes'][i].append([nodes[-1], nodes[0]]) 69 else: # the last is the odd one out 70 rule['nodes'][i].append([nodes[0], nodes[-1]]) 71 else: 72 X = [x, y, mpfr(1)] 73 B = biunit_to_bary(X) 74 rule['nodes'][i].append(B) 75 rules[power] = rule 76 rules["type"] = "triangle" 77 rules["prefix"] = "WVTri" 78 rules["url"] = "https://doi.org/10.1016/j.camwa.2015.03.017" 79 return rules 80 81 82def tetrahedron_rules(rules_path): 83 """ 84 Gather and format quadrature rule data from the unpacked zip file found 85 at https://doi.org/10.1002/nme.6528 86 """ 87 88 rules = {} 89 rules_path = os.path.join(rules_path, "cubatures_tet_sym_repository") 90 for file in glob.glob(os.path.join(rules_path, "*_compact.txt")): 91 name = os.path.basename(file) 92 descrip = os.path.splitext(name)[0] 93 parts = re.search('cubature_tet_sym_p([0-9]+)_n([0-9]+)_compact', descrip) 94 power = int(parts.group(1)) 95 num_points = int(parts.group(2)) 96 counts = () 97 with open(os.path.join(rules_path,file),"r") as f: 98 count_regex = re.findall('=([0-9]+)', f.readline()) 99 counts = tuple(int(s) for s in count_regex) 100 rule = {} 101 rule['num_nodes'] = num_points 102 rule['counts'] = counts 103 rule['nodes'] = [[], [], [], [], []] 104 rule['weights'] = [] 105 exname = name.replace('compact','expand_baryc') 106 with open(os.path.join(rules_path, exname),"r") as f: 107 for (i,(count, multiplicity)) in enumerate(zip(counts, (1,4,6,12,24))): 108 for c in range(count): 109 for m in range(multiplicity): 110 w, x, y, z, weight = tuple(mpfr(s) for s in re.split('\s+',f.readline().strip())) 111 if m == 0: 112 rule['weights'].append((weight * 4) / 3) 113 if i == 0: # (a, a, a, a) 114 rule['nodes'][i].append([w]) 115 elif i == 1: # (a, a, a, b) 116 nodes = sorted([w, x, y, z]) 117 diffs = [b - a for (a,b) in zip(nodes[:-1],nodes[1:])] 118 if diffs[0] > diffs[-1]: # first index is the odd one out 119 rule['nodes'][i].append([nodes[-1], nodes[0]]) 120 else: # last index is the odd one out 121 rule['nodes'][i].append([nodes[0], nodes[-1]]) 122 elif i == 2: # (a, a, b, b) 123 nodes = sorted([w, x, y, z]) 124 rule['nodes'][i].append([nodes[0], nodes[-1]]) 125 elif i == 3: # (a, a, b, c) 126 nodes = sorted([w, x, y, z]) 127 diffs = [b - a for (a,b) in zip(nodes[:-1],nodes[1:])] 128 j = 0 129 for k in range(len(diffs)): 130 if diffs[k] < diffs[j]: 131 j = k 132 if j == 0: # 0 and 1 are the same 133 rule['nodes'][i].append([nodes[0], nodes[2], nodes[3]]) 134 elif j == 1: # 1 and 2 are the same 135 rule['nodes'][i].append([nodes[1], nodes[0], nodes[3]]) 136 else: # 2 and 3 are the same 137 rule['nodes'][i].append([nodes[2], nodes[0], nodes[1]]) 138 else: 139 rule['nodes'][i].append([w, x, y, z]) 140 rules[power] = rule 141 rules["type"] = "tetrahedron" 142 rules["prefix"] = "JSTet" 143 rules["url"] = "https://doi.org/10.1002/nme.6528" 144 return rules 145 146 147def create_header(rules): 148 powers = [i for i in filter(lambda x: type(x) == int, rules.keys())] 149 max_power = max(powers) 150 filename = f"petscdt{rules['type'][:3]}quadrules.h" 151 guard = filename.replace('.','_').upper() 152 lines = [] 153 lines.append(f"""// Minimal symmetric quadrature rules for a {rules['type']} from {rules['url']} 154// This file was generated by lib/petsc/bin/maint/petscdt_create_quadrature_headers.py 155#if !defined({guard}) 156#define {guard} 157 158#include <petscsys.h> 159""") 160 printer_wide = pprint.PrettyPrinter(width=260) 161 printer_narrow = pprint.PrettyPrinter(width=58) 162 prefix = f"PetscDT{rules['prefix']}Quad" 163 for key in sorted(powers): 164 rule = rules[key] 165 lines.append(f"static const PetscReal {prefix}_{key}_weights[] = {{"); 166 weights = printer_narrow.pformat(["PetscRealConstant({0:e})".format(w) for w in rule['weights']]).replace('[',' ').replace("'","").strip(']') 167 lines.append(f"{weights}") 168 lines.append("};") 169 lines.append("") 170 lines.append(f"static const PetscReal {prefix}_{key}_orbits[] = {{"); 171 for (i,nodes) in enumerate(rule['nodes']): 172 for (j,node) in enumerate(nodes): 173 fmt = printer_wide.pformat(["PetscRealConstant({0:e})".format(w) for w in node]).replace('[',' ').replace("'","") 174 if sum([len(s) for s in rule['nodes'][i+1:]]) == 0 and j == len(nodes) - 1: 175 fmt = re.sub('\s*]','',fmt) 176 else: 177 fmt = re.sub('\s*]',',',fmt) 178 lines.append(fmt) 179 lines.append("};") 180 lines.append("") 181 182 lines.append(f"static const PetscInt {prefix}_max_degree = {max_power};") 183 lines.append("") 184 degree_to_power = [-1 for i in range(max_power+1)] 185 for i in range(max_power, -1, -1): 186 if i in rules.keys(): 187 degree_to_power[i] = i 188 else: 189 degree_to_power[i] = degree_to_power[i+1] 190 191 num_nodes = [rules[degree_to_power[d]]['num_nodes'] for d in range(max_power+1)] 192 lines.append(f"static const PetscInt {prefix}_num_nodes[] = {{") 193 lines.append(f"{printer_narrow.pformat(num_nodes).replace('[',' ').strip(']')}") 194 lines.append("};") 195 lines.append("") 196 197 num_orbits = [rules[degree_to_power[d]]['counts'] for d in range(max_power+1)] 198 lines.append(f"static const PetscInt {prefix}_num_orbits[] = {{") 199 lines.append(f"{printer_narrow.pformat(num_orbits).replace('[',' ').strip(']').replace('(','').replace(')','')}") 200 lines.append("};") 201 lines.append("") 202 203 for suff in ["weights", "orbits"]: 204 lines.append(f"static const PetscReal *{prefix}_{suff}[] = {{"); 205 for i in range(max_power+1): 206 if i == max_power: 207 lines.append(f" {prefix}_{degree_to_power[i]}_{suff}"); 208 else: 209 lines.append(f" {prefix}_{degree_to_power[i]}_{suff},"); 210 lines.append("};") 211 lines.append("") 212 213 lines.append(f"#endif // #define {guard}") 214 215 with open(filename, "w") as f: 216 f.writelines('\n'.join(lines) + '\n') 217 218 219if __name__ == '__main__': 220 import tempfile 221 import requests 222 from zipfile import ZipFile 223 from tarfile import TarFile 224 import argparse 225 import pathlib 226 227 parser = argparse.ArgumentParser(description="Create quadrature headers from published compact representations") 228 parser.add_argument('--tri', type=pathlib.Path) 229 parser.add_argument('--tet', type=pathlib.Path) 230 args = parser.parse_args() 231 232 def unzip_and_create_header(file, extracter, rules): 233 with tempfile.TemporaryDirectory() as tmpdirname: 234 _dir = os.path.join(tmpdirname, "rules") 235 with extracter(file) as _archive: 236 _archive.extractall(path=_dir) 237 create_header(rules(_dir)) 238 239 if args.tri: 240 unzip_and_create_header(args.tri, ZipFile, triangle_rules) 241 242 if args.tet: 243 unzip_and_create_header(args.tet, TarFile, tetrahedron_rules) 244 245